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1.0  SUMMARY 


This  work  models  charging  of  space  objects  and  investigates  the  effect  it  has  on  the  orbital  evo¬ 
lution  of  those  object.  In  the  approach  that  is  presented  in  this  research,  a  spherical  conductor 
model  has  been  adopted.  In  a  steady  state  approach,  charging  levels  have  been  derived  Under 
Low-charging  plasma  condition  and  high-charging  plasma  condition.  A  dipole  model  for  the  Earth 
magnetosphere  has  been  adopted.  Formulations  for  Luni-Solar  gravitational  force  and  solar  radia¬ 
tion  pressure  are  taken  from  Frueh  et  al.  [1],  In  this  research,  we  include  Earth  gravitational  spher¬ 
ical  expansion  up  to  order  and  degree  12  [2].  Gosudarstvennyy  Standart  (GOST)  model  [3]  is  used 
to  simulate  atmospheric  drag  force  in  low  Earth  orbit.  Atmospheric  drag  is  assumed  to  be  absent  in 
geosynchronous  orbits.  This  research  will  improve  understanding  of  space  plasma  environment  in 
the  low  Earth  orbits  (LEO)  and  geosynchronous  Earth  orbits  (GEO),  and  how  the  plasma  interacts 
with  space  objects  leading  to  charge  development.  It  will  also  shed  light  into  photo-emission  effect 
on  charging.  The  research  will  investigate  effects  of  Lorentz  force  on  orbital  evolution  of  standard 
low  area  to  mass  ratio  (LAMR)  objects  and  highly  perturbation  susceptible  high  area  to  mass  ratio 
(HAMR)  objects.  Furthermore,  it  will  shed  light  into  effects  of  presence/absence  of  sunlight  and 
orbit  inclination  on  the  evolution  of  orbital  parameters. 

2.0  INTRODUCTION 

High  fidelity  propagation  models  are  integral  for  the  long  and  short  term  orbital  prediction  of  near 
Earth  space  objects.  For  a  robust  space  catalog,  the  orbit  prediction  of  an  individual  object  has  to  be 
better  than  half  the  field  of  view  of  the  observing  sensor;  which  is  of  the  order  of  several  degrees 
down  to  less  than  half  a  degree.  The  frequency  of  observations  can  be  really  low,  sometimes 
the  time  between  two  consecutive  observations  can  exceed  30  orbital  periods.  Ideally,  orbital 
predictions  need  to  be  on  the  arc  second  level  several  days  ahead  to  allow  for  an  informed  decision 
about  possible  collision  avoidance  maneuvers. 

In  the  orbit  prediction,  usually  perturbation  forces  include  higher  harmonics  of  Earth’s  gravity, 
Luni-Solar  gravitational  forces,  solar  radiation  pressure  (cannonball  model),  and  atmospheric  drag. 
An  important  question,  however,  remains  is  whether  the  inclusion  of  other  perturbation  forces 
beyond  those  mentioned  above  is  significant  for  precise  orbit  propagation.  In  this  research,  Lorentz 
force  has  been  considered  as  an  additional  perturbation. 

A  body  obtains  a  charge  when  there  is  a  difference  between  the  fluxes  of  electrons  and  ions  received 
by  the  body.  In  this  case  the  object  develops  a  potential  so  as  to  reduce  the  flux  of  more  available 
particle  and  increase  the  flux  of  less  available  particle  until  both  of  them  equals  out  and  the  body 
reaches  an  equilibrium  state.  Natural  space  plasma  is  sufficient  to  introduce  significant  charges. 
Spacecraft  charging  research  has  first  come  into  focus  because  of  several  satellite  anomalies  ob¬ 
served  during  the  1970s.  The  1973  failure  of  Defense  Satellite  Communications  System  (DSCS) 
9431  satellite  of  US  Air  Force  because  of  a  discharge  lead  National  Aeronautics  and  Space  Admin¬ 
istration  (NASA)  and  the  US  Air  Force  to  investigate  the  effects  of  charging  and  develop 
technologies  to  mitigate  the  same.  Accumulated  charges  may  result  in  currents  entering  the 
tender  circuits  and  impeding  their  functionality.  Sometimes,  differential  charging  across  various 
parts  of  the  space-craft  can  be  high  enough  to  initiate  arcing.  The  effects  of  arcing  on  the  on¬ 
board  electronics  can 
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vary  from  mild  to  fatal  depending  on  the  magnitude  of  differential  voltage  and  the  spacecraft  ma¬ 
terial.  Charged  particles  can  also  create  strong  electric  fields,  which  can  interfere  with  sensitive 
scientific  instruments,  resulting  in  wrong  data  measurements. 

Several  studies  investigating  charging  in  the  near  geosynchronous  orbit  have  been  done.  DeForest 
[4]  investigated  charging  of  ATS-5  satellite  in  geosynchronous  orbit.  He  found  that  there  were 
high  levels  of  charging  of  the  order  of  thousands  of  negative  volts  during  eclipse  and  few  hundred 
negative  volts  during  sunlit  conditions.  The  difference  in  charging  levels  between  shadow  and  sun 
light  conditions  result  from  the  emission  of  photo  electrons.  Mulen  et  al.  [5]  analyzed  the  high- 
level  charging  under  sunlit  condition.  Their  work  was  based  on  data  from  Spacecraft  Charging 
at  High  Altitudes  (SCATHA)  satellite.  Although  less  common,  high-level  charging  can  also  be 
experienced  in  low  Earth  orbit.  Defense  Meteorological  Satellite  Program  (DMSP)  spacecraft  in 
low  Earth  polar  orbit  was  found  to  charge  up  to  voltage  as  high  as  -2000  volts  [6].  In  addition  to 
the  few  mentioned  studies,  comprehensive  texts  by  Shu  [7],  Hastings  and  Garrett  [8]  have  greatly 
furthered  theoretical  knowledge  of  space  environment  and  charging. 

Several  existing  software  are  available  for  simulating  spacecraft  charging.  Jointly  developed  by 
NASA  and  Air  Force  Research  Laboratory  (AFRL),  NASA/Air  Force  Spacecraft  Charging  Anal¬ 
ysis  Program  (NASCAP-2K)  [9]  is  probably  the  most  widely  known  software  in  this  regard.  The 
computations  employed  by  NASCAP-2K  involves  both  analytic  and  particle-in-cell  (PIC)  meth¬ 
ods,  and  it  is  based  on  previously  developed  NASA  Charging  Analyzer  Program  for  Low-Earth 
Orbit  (NASCAP-LEO),  NASA  Charging  Analyzer  Program  for  Geosynchronous  Orbit  (NASCAP- 
GEO)  and  Potentials  Of  Large  objects  in  the  Auroral  Region  (POLAR).  The  Spacecraft  Plasma 
Interaction  System  (SPIS)  is  an  open-source  software  under  development  and  it  is  based  on  three- 
dimensional  particle-in-cell  computations.  It  includes  several  detailed  interactions  like  photoelec¬ 
tric  effect,  secondary  emission,  etc.  but  the  software  graphics  and  documentation  still  need  to  be 
improved  so  as  to  increase  user-friendliness.  Commercially  available  spacecraft-charge  simulating 
software  include  Multi-Utility  Spacecraft  Charging  Analysis  Tool  (MUSCAT)  [10]  and  Spacecraft 
Charging  Software  (SPARCS)  [11]. 

The  focus  of  the  aforementioned  tools  is  the  precise  investigation  of  charging  levels  with  respect 
to  a  pre-launch  satellite  design.  The  interest  is  to  compute  charge  spikes  for  a  given  orbital  design. 
Hence,  the  calculations  are  based  upon  the  assumption  that  the  objects  are  placed  on  Keplerian 
orbits. 

However,  when  the  charged  object  moves  relative  to  the  Earth  magnetosphere,  Lorentz  forces 
are  induced.  Those  forces  result  in  orbital  perturbations.  First  steps  towards  evaluating  how  signif¬ 
icant  Lorentz  force  perturbations  are  in  the  orbital  evolution  of  natural  objects,  a  few  preliminary 
studies  have  been  made.  Antal  et  al.  [12]  studied  the  dynamics  of  charged  micron  and  sub-micron 
sized  space  debris  particles  in  Earth’s  plasma  environment.  Aziz  [13]  used  Lagrange  planetary 
equations  to  derive  analytical  formulas  for  perturbations  in  orbital  elements  of  middle  Earth  LA- 
GEOS  satellite  due  to  Lorentz  force,  while  assuming  a  fixed  value,  rather  than  the  one  induced  by 
the  actual  plasma  conditions  in  the  orbital  region,  for  body  charge.  Peng  et  al  [14]  derived  Gauss 
variational  equations  for  Lorentz  force  perturbations  under  the  assumption  of  an  invariant  Earth 
magnetic  dipole.  Frueh  et  al.  [15]  simulated  orbital  evolution  of  HAMR  flat  plates  in  geostation¬ 
ary  orbit.  They  took  into  account  solar  radiation  pressure,  Earth’s  gravity,  third  body  perturbations 
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and  Lorentz  force.  In  their  simulation,  charging  plasma  values  were  computed  in  interaction  of 
the  orbital  propagation  and  recurrent  use  of  the  NASCAP-2k  tool,  to  allow  the  computation  of  the 
charging  in  the  natural  and  changing  plasma  environment. 

In  this  study,  the  work  from  Frueh  et  al  has  been  developed  further.  Instead  of  relying  on  the 
NASCAP  tool,  the  natural  charging  effects  have  been  computed  alongside  the  orbital  propagation, 
in  order  to  be  able  to  investigate  the  effect  of  the  Lorentz  perturbations  as  realistically  as  possible 
in  a  computationally  feasible  integrated  fashion.  In  order  to  allow  for  computational  fast  computa¬ 
tions,  the  analytic  variational  equations  form  Peng  et  al  [14]  have  been  expanded. 

This  report  is  organized  as  the  following:  In  the  next  section,  the  model  and  assumptions  are 
introduced.  The  model  consist  of  the  dipole  model  of  the  Earth  magnetosphere  with  time  varying 
coefficients,  the  plasma  environment.  Subsequently,  the  modeling  of  the  space  object  charging 
and  its  capacitance  is  shown.  The  analytic  variational  equation  expressions  are  then  introduced. 
In  the  next  section  the  numerical  test  cases  are  shown  and  discussed.  Test  orbits  of  high  and  low 
area-to-mass  ratio  objects  in  low  Earth  orbit  and  geosynchronous  orbital  region  have  been  propa¬ 
gated.  The  report  concludes  with  a  summary  of  the  findings  and  the  conclusions.  In  Appendices 
A  and  B,  coefficients  of  the  variational  equations  developed  in  this  research  can  be  found,  and  in 
Appendix  C,  Peng’s  [14]  variational  equations  are  listed  for  comparison.  Symbols,  abbreviations 
and  acronyms  follow  the  appendices. 


3.0  METHODS,  ASSUMPTIONS  AND  PROCEDURES 
3.1  Assumptions 

In  the  approach  that  is  presented  in  this  research,  a  spherical  conductor  model  has  been  adopted.  In 
a  steady  state  approach,  charging  levels  have  been  derived  Under  Low-charging  plasma  condition 
and  high-charging  plasma  condition.  A  dipole  model  for  the  Earth  magnetosphere  with  time  vary¬ 
ing  coefficients  has  been  adopted  [  16].  Formulations  for  Luni-Solar  gravitational  force  and  solar 
radiation  pressure  are  taken  from  Frueh  et  al.  [1],  Earth  gravitational  spherical  expansion  up  to 
order  and  degree  12  [2]  are  modeled.  Model  [3]  is  used  to  simulate  atmospheric  drag  force  in 
low  Earth  orbit.  Atmospheric  drag  is  assumed  to  be  absent  in  geosynchronous  orbits. 


3.2  Orbital  Differential  Equation 

A  three  degree  of  freedom  motion  is  assumed  for  the  debris  since  a  spherical  geometry  is  consid¬ 
ered  in  this  study.  The  following  equation  of  motion  holds  for  center  of  mass, 


x  =  —GMqVV (x)  —  G  £  Mk 

k=  1,2 


x-xk  Xk 
x-xk\3  xl 


+  52^/ 


(1) 


where  x  is  the  geocentric  position  of  the  object,  G  the  gravitational  constant,  M®  the  Earth  mass  and 
V  (x)  the  Earth  gravitational  potential.  The  formulation  of  Pines  [2]  has  been  implemented  for  the 
gravitational  potential  representations,  and  the  user  can  choose  order  and  degree  representations 
up  to  100.  The  third  body  gravitational  perturbations  of  the  Sun  and  Moon  (k=l,  2)  with  the 
states  xk  are  modeled  via  their  center  of  mass.  The  Sun  and  Moon  masses  are  represented  by  Mk 


3 

Approved  for  public  release;  distribution  is  unlimited. 


(k=l,  2).  Finally,  J^cq  is  the  sum  over  all  non- gravitational  accelerations  acting  on  the  satellite. 
Solar  radiation  pressure,  atmospheric  drag  and  Lorentz  force  are  the  three  non-gravitational  forces 
considered  in  this  paper. 


The  acceleration  due  to  solar  radiation  pressure  for  a  spherical  body  of  radius  r  is  modeled  as, 


amd  — 


47 rr2  E  A2  / 1  1  \  * 

- F - (  A  Q^d  )  ^ 

m  c  \x  —  Xq\z  \4  9  / 


(2) 


where  m  represents  the  total  mass  of  the  body,  E  represents  solar  flux  at  1  AU,  c  represents  speed 
of  light,  A©  represents  astronomical  unit  and  Q  is  the  diffuse  reflection  coefficient,  x,  xQ,  S 
represent  geocentric  position  of  object,  geocentric  position  of  sun  and  direction  of  radiation  source, 
respectively. 


The  acceleration  due  to  atmospheric  drag  force  is  modeled  as, 

_  _  _  1  cdA  re  y P  |  Vre/ v  |  V rei v 

®drae  ~  X  wJ 

2  m 

where  cd  is  drag  coefficient,  Aref  is  reference  area,  p  is  atmospheric  density  and  v,e/v  is  velocity 
of  body  relative  to  atmosphere.  Aref  is  taken  as  spherical  cross-sectional  area  i.e.,  nr2  and  the 
atmosphere  is  assumed  to  be  co-rotating  with  Earth.  The  density  p  is  obtained  using  GOST  model 
as, 


p  =  pnkokik2k3k4  (4) 

where  the  parameters  p„  and  kj  (j=  1,2, 3, 4)  are  semi-empirical  factors  that  take  into  account  various 
aspects  like  diurnal  variation,  solar  cycle,  semi-annual  variation,  etc.  and  can  be  obtained  from 
Vallado  [3], 


The  acceleration  due  to  Lorentz  force  is  modeled  as, 


a charge 


qVrel  X  B 


m 


(5) 


where  B  is  the  Earth  magnetic  field,  q  is  the  surface  charge  of  spherical  body  and  vrei  is  the  velocity 
of  the  object  relative  to  the  magnetic  field.  Charge  q  can  be  obtained  from  potential  <j>  of  the  space 
object  using  the  following  relation, 

q  =  C$  (6) 

where  C  is  body  capacitance.  The  body  potential  (j)  is  obtained  by  solving  the  so-called  current 
equation. 


A  variety  of  currents  contribute  to  body  charging.  A  body  in  space  acts  like  a  node  in  an  electrical 
circuit  and  hence  Kirchhoff’s  current  law  is  applicable.  At  equilibrium,  according  to  Kirchhoff’s 
law, 

LW)  =  0  (7) 

j 

li  3“  le  T"  Iph  T  Ce  3“  Ibsc  T  hnisc  —  0  (8) 
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where  /,  is  plasma  ion  current,  Ie  is  plasma  electron  current,  Iph  is  photoelectric  current,  Ise  is 
secondary  electron  current,  4SC  is  backscattered  current  and  /m(SC  is  miscellaneous  current.  Mis¬ 
cellaneous  current  includes  all  currents  apart  from  ones  just  mentioned,  and  it  will  be  neglected 
in  all  the  theoretical  developments  and  simulations  in  this  work.  Some  important  components 
constituting  miscellaneous  current  may  include  electron  beam  emission,  neutral  particle  generated 
secondary  electrons  and  ions,  etc.  0  appearing  in  Equation  (7)  is  body  potential. 


In  non-equilibrium  scenario,  charging  takes  place  according  to  the  differential  equation,  Y*jlj  = 
C^.  The  time  scale  for  charging  of  a  spherical  conductor  can  be  given  approximately  by  [25], 


C</> 

Anr2J 


(9) 


where  the  subscript  ts  stands  for  transient  time.  The  parameters  C,  0,  r,  J  represent  body  capac¬ 
itance,  body  potential,  body  radius  and  ambient  flux,  respectively.  In  this  research,  the  charging 
time  scale  is  much  smaller  compared  to  orbit  propagation  time  steps  and  hence  equilibrium  current 
condition  is  assumed. 


3.3  Earth  Magnetosphere 

Lorentz  acceleration  given  by  Equation  (5)  depends  upon  Earth’s  magnetosphere.  Typically, 
magneto-sphere  models  can  be  grouped  under  two  classes:  statistical  and  physics-based.  The 
Physics-based  models  are  computationally  expensive  and  difficult  to  solve  as  they  require 
solving  numerically  several  nonlinear  partial  differential  equations  on  a  three-dimensional  grid 
of  points.  They  take  magnetohydrodynamic  (MHD)  flows  of  the  plasma  into  consideration  for 
simulating  the  magnetic  field.  Statistical  models,  on  the  other  hand,  are  based  on  empirical 
formulations  and  are  most  reliable  under  weak  storm  scenarios. 

Early  geomagnetosphere  models  were  primarily  physics-based  due  to  unavailability  of  sufficient 
space-data,  but  over  the  last  three  decades,  with  the  availability  of  huge  amount  of  data  from 
satellites  and  space-instruments,  trend  has  shifted  towards  use  of  statistical  models.  Mead-Fairfield 
geomagnetic  field  model  [17]  is  one  of  the  earliest  statistical  models.  The  model  is  valid  till  17 
earth  radii,  and  it  is  created  with  the  help  of  data  from  Interplanetary  Monitoring  Platform  (IMP) 
satellites  spanning  over  the  period  from  1966  to  1972.  The  Mead-Fairfield  model  is  linear  in  tilt, 
quadratic  in  position  and  power  series  (second  order)  in  solar  magnetic  coordinates.  Following 
this,  Tsyganenko,  Usmanov,  and  Sitnov  [18,  19]  contributed  significantly  to  the  development  of 
statistical  models.  The  Tsyganenko  magnetic  field  model  is  based  on  data  from  various  satellites 
like  IMP,  Highly  Eccentric  Orbit  Satellite  (HEOS),  International  Sun-Earth  Explorer  (ISEE),  Polar, 
Geotail,  etc.  The  model  takes  into  consideration  the  effect  of  ring  current,  magnetotail  current 
system,  magnetopause  current,  field-aligned  current  system. 

Olson-Pfitzer  q  uiet,  1  977  [  20]  i  s  a  nother  e  arly  s  tatistical  E  arth  m  agnetic  fi  eld  mo  del  developed 
with  the  aid  of  data  from  Orbiting  Geophysical  Observatory  (OGO)  3  and  5.  The  model’s  va¬ 
lidity  extends  from  dayside  subsolar  magnetosphere  to  farther  than  lunar  orbit  in  the  nightside 
magnetotail.  This  model,  however,  applies  to  only  quiet  conditions.  It  includes  the  contributions 
from  tail,  magnetopause  and  ring  currents  and  the  core  field  is  approximated  by  a  fixed  dipole. 
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A  modification  of  Olson-Pfitzer  Field  Model  (1977)  is  Olson-Pfitzer  Dynamic  Model,  1988 

[21] .  It  takes  into  consideration  the  variation  in  the  size  and  strength  of  primary 
magnetospheric  current  systems  resulting  from  interplanetary  sources.  The  model  is 
appropriate  for  dayside  magnetic  field  at  geosynchronous  orbit  for  any  magnetic  condition 
and  nightside  field  for  quiet  conditions.  The  Ostapenko-Maltsev  Magnetic  Field  Model,  1997 

[22]  is  another  model  which  uses  the  database  of  Fairfield  et  al.  [23]  and  valid  for  3 Re  to 
IRe  .  The  model  is  dependent  on  geomagnetic  indices  Kp  and  Dst  ,  solar  wind  dynamic 
pressure  and  z-component  of  interplanetary  magnetic  field. 


The  Lyon-Fedder-Mobarry  (LFM)  code  is  a  three-dimensional  MHD  modeling  of  the  magneto¬ 
sphere  and  the  work  on  the  code  first  started  around  1985.  Over  the  years,  the  code  has  been 
continuously  improved  to  accommodate  better  resolving  power  transport;  an  adaptable  grid  based 
on  the  problem  statement;  addressing  the  problem  of  non-zero  V  •  B:  improvement  upon  the  cases 
of  high  Alfven  speed,  low  ratio  of  plasma  pressure  to  magnetic  pressure  and  strong  magnetic  field 
gradients;  and  integration  of  ionospheric  model  [24]. 

The  work  carried  out  here  uses  the  empirical  dipole  magnetic  field  model.  The  model  is  a  first 
order  approximation  of  Earth’s  actual  magnetic  field.  The  implementation  of  dipole  model  is 
simple,  and  is  particularly  accurate  at  low  Earth  altitudes.  Equations  (10),  (11),  (12)  represent  the 
dipole  model  [16],  where  the  parameters  17,  £  represent  geocentric  distance,  co-elevation  and 

East  longitude  from  Greenwich,  respectively,  /y,  is  equatorial  radius  of  Earth.  The  parameters 
gj  and  hj  are  International  Geomagnetic  Reference  Field  (IGRF)  coefficients  and  can  be  obtained 
from  data  released  by  International  Association  of  Geomagnetism  and  Aeronomy  (IAGA). 


Bn  =2 

Bs  =  ( 


V 
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3.4  Space-Plasma  Environment 

Body  charge  depends  on  longitude,  local  time  and  altitude,  because  of  varying  plasma  conditions 
[7].  Current  work  focuses  on  low  Earth  altitudes  (approximately  150  Km  to  2000  Km)  and  geosyn¬ 
chronous  altitude  only.  In  the  current  work,  for  the  low  Earth  regions,  local  time  and  longitudinal 
variations  of  plasma  data  are  neglected  focusing  only  on  altitude  dependent  effects  as  they  ap¬ 
pear  to  be  the  most  crucial  ones.  For  the  geosynchronous  region,  only  the  local  time  variations 
are  considered  because  of  limitations  in  data  availability.  Plasma  at  low  Earth  orbits  is  usually  of 
low  energy  (typically  less  than  .1  eV  energy)  and  high  density  (typically  105  cm-3  or  higher).  At 
geosynchronous  altitudes,  plasma  density  is  relatively  small  (around  1  cm~ 3),  but  plasma  energy  is 
high,  typically  around  100  eV,  often  reaching  KeV  ranges  during  geomagnetic  storm  activities.  A 
Geomagnetic  storm  is  a  space  weather  phenomena  resulting  from  rapid  increase  in  energy  transfer 
from  Sun  to  Earth’s  magnetosphere  and  typically  caused  by  solar  coronal  mass  ejections  (CME). 
Bodies  in  geosynchronous  orbit  typically  get  positively  charged  in  sunlight  due  to  emission  of 
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photo-electrons.  The  photoelectric  current  in  geosynchronous  region  is  of  similar  order  as  that 
of  the  dominant  plasma  electron  current  and  hence  it  is  a  major  contributor  to  body  charge.  For 
low  Earth  region,  photoelectric  current  is  relatively  insignificant  as  plasma  currents  are  much  larger 
than  photoelectric  current.  In  geosynchronous  shadow  region  and  in  all  of  low  Earth  region,  bodies 
are  typically  negatively  charged  because  of  dominant  plasma  electron  current. 

Once  the  voltage  is  determined  from  current  equilibrium,  the  charge  can  be  determined  using  ca¬ 
pacitance.  The  capacitance  is  both  object  geometry  dependent  and  plasma  dependent.  The  plasma 
dependence  is  characterized  by  means  of  the  so-called  Debye  length.  When  a  body  gets  charged, 
it  attracts  the  species  of  opposite  polarity  and  repels  the  species  of  same  polarity  and  as  a  result, 
a  layer  adjacent  to  body  is  formed  which  is  rich  in  opposite  polar  species  and  poor  in  same  polar 
species.  A  Debye  length  is  roughly  the  distance  outside  which  no  effect  of  the  body  charge  is  felt 
i.e.,  the  electric  field  of  the  body  roughly  cancels  out  the  electric  field  of  the  opposite  polar  species 
surrounding  the  body.  Since  plasma  density  is  less  in  geosynchronous  Earth  orbit  as  compared  to 
low  Earth  orbit,  more  thickness  of  the  opposite  polar  species  would  be  required  in  geosynchronous 
orbit  to  cancel  out  the  electric  field  due  to  body  charging  and  hence  Debye  length  is  comparatively 
large  or  sheath  is  thick.  Opposite  is  true  for  low  Earth  orbit  and  the  plasma  sheath  is  thin.  For 
geosynchronous  orbit,  Debye  length  is  typically  greater  than  200m  whereas  for  low  Earth  orbit,  it 
is  in  millimeter  range.  As  will  be  seen  later,  formulation  of  current  depends  on  sheath  thickness. 
It  is  worth  mentioning  that  in  terms  of  analytical  convenience,  sphere,  infinite  plate  and  infinite 
cylinder  are  three  geometries  which  are  usually  looked  into,  because  of  the  symmetry  in  their  field. 

Plasma  electron  and  ion  currents,  secondary  electron  and  ion  currents,  backscattered  electron  cur¬ 
rent  and  photoelectric  current  are  modeled  using  steady  state  approach. 

3.5  Plasma  Current  Model 

Depending  on  whether  a  particular  kind  (polarity)  of  plasma  species  is  attracted  or  repelled,  plasma 
current  expression  can  be  different. 

This  section  derives  analytic  expression  [26],  [27]  for  flux  of  particles  repelled  by  the  polarity 
of  the  space  object  in  low  Earth  and  geosynchronous  orbits. 

Let  us  consider  a  ID  system  of  particles  as  shown  in  Figure  1.  Let  n  be  the  number  of  particles 
per  unit  length.  For  any  particular  speed  v  and  a  time  span  t,  only  the  particles  within  a 
distance  of  vt  from  the  wall  will  be  able  to  reach  it.  Thus,  if  v~  represents  the  mean  speed,  then 
the  average  number  of  particles  reaching  the  wall  in  a  time  span  of  t  is  given  as, 

A  inear  density ,  n 

Nr  —  ( strike  radius )  ( - - - -J  —  vt  ■  -  (13) 

And  the  particle  strike  rate  is  given  as, 

Nr  vn 

T  =  Y 
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Half  Moving  Right 
and  the  Other  ^ 
Half  Moving  Left 


Figure  1.  ID  System  of  Particles  Hitting  a  Wall  [27] 

This  idea  can  be  extended  to  3D.  Let  f(v)  be  the  probability  density  function  for  velocity  of  par¬ 
ticles,  then  the  mean  number  of  particles  per  unit  volume  between  velocities  v  and  v  +  dv  is  given 
by  f(v) d3v  =  f(v)  sin  0d0d\j/v2dv,  where  v  is  velocity  magnitude  and  0,  yr  represent  polar  coor¬ 
dinates. 

With  reference  to  Figure  2,  number  of  particles  having  velocity  ~v  reaching  the  elemental  area  A 
in  time  span  t  can  be  written  as, 

Nr  —  (No.  per  unit  volume)  (volume  of  slanted  cylinder)  —  [f(v)d3v]  [Avtcos  0]  (14) 


vt 


Area  A 


Figure  2.  Side  View  of  Wall  Area  A  Being  Hit  by  Particles  [27] 

Number  of  particles  hitting  wall  per  unit  area  per  unit  time  with  velocity  v  is, 

—  =  f(v)d3vvc  os  0 
At  w 

To  obtain  the  flux  0o,  the  above  expression  needs  to  be  integrated  over  all  velocity  ranges, 


(j) o  —  J  f(v)v  cosOd3  v 

(15) 

(f>o—  j  j  j  f(v)v cos 9 sm9dddyrv2dv 

(16) 

C  1  n  ^  r27l 

(j)o—  Jv3f(v)dv--J  sin2 9d9  ■  j  d\jf 

(17) 

(j)o  =  n  1  v3f(v)dv 

(18) 
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If  the  particles  follow  Maxwell-Boltzmann  distribution,  then  the  plasma  velocity  distribution,  f(v) 
and  the  plasma  speed  distribution,  F(v)  are  given  as, 


/(v)  =  n 


2nkT 


e  2 kT 


F{y)  =  Ann 


2nkT 


3 

2  2 
_ mv 

v  e 


where  n,  m,  k,  T  represent  particle  number  density,  particle  mass,  Boltzmann  constant  and  plasma 
distribution  temperature,  respectively. 


Using  Equations  (19),  (20),  the  flux  in  Equation  (18)  can  then  be  written  in  terms  of  distribution 
function  as, 

o  =  -7--~  /  vF(v)dv  (21) 

An  n  .1 


,  n  ZkT  /  kT 

~  4  V  Itiii  iKni  <23) 

where  v  is  mean  speed.  When  the  spacecraft  is  charged  to  potential  V,  the  probability  distribution 
changes  from  f(E)  to  f(E  +  qV)  [7],  where  q  is  particle  charge  and  E  is  particle  energy.  Thus,  the 
flux  is  given  as, 


m  \  2  i"(v2+4f)  ,  1 


2  nkT 


e  2*  t  v2 -sin26  dOdydv 


qV 

(j)  —  e  k  i  (j)o 

where  (j) o  is  the  flux  obtained  earlier  in  Equation  (23)  in  absence  of  spacecraft  charging. 
The  current  density  to  an  uncharged  body  can  be  written  in  terms  of  flux  as, 

,  /  kT  1  fl kT 

and  the  current  density  to  a  body  charged  to  a  potential  V  is  given  as, 


qV  qV 

J  =  q(j)  —  e  kr  q(j)Q  =  e  kr  JQ 


It  is  to  be  noted  again  that  the  above  derivation  is  based  on  the  assumption  that  the  plasma  particles 
follow  Maxwell-Boltzmann  distribution  under  thermal  equilibrium  condition.  The  assumption  is 
usually  valid  for  calm  environment.  In  actuality,  no  single  distribution  can  accurately  model  the 
plasma  particles  because  of  various  disturbances. 
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The  flux  for  particles  can  be  discriminated  by  the  attracted  by  debris  polarity.  In  the  geosyn¬ 
chronous  region,  where  the  plasma  density  is  low,  the  incoming  current  density  collected  by  a 
perfectly  absorbing  spherical  or  infinitely  long  cylindrical  body  from  a  collisionless,  stationary, 
isotropic  plasma  is  limited  by  the  particles’  orbital  angular  momenta.  This  current  is  known  as 
orbit-limited  current  [28].  For  orbit-limited  sphere,  for  an  attracted  particle  approaching  the  body 
from  infinity,  conservation  of  energy  and  angular  momentum  results  in, 

\mvl  =  ^wv(rs)2  +  qV  (rs)  (28) 

mR\v  o  =  mrsv(rs)  (29) 

where  vo  is  velocity  in  ambient  medium,  v  is  velocity  at  object  surface,  R\  is  impact  parameter,  V 
is  body  potential  and  rs  is  object  radius.  Only  particles  having  radius  less  than  R\  will  reach  rs. 
Manipulations  of  Equations  (28)  and  (29)  would  produce, 

1  2  1  [RlVo\2  2q 

-mv0  =  -m^— j  +-V(rs)  (30) 


R2  =  r2 
rs 


1- 


2  qV  (rs 


mvn 


(31) 


(R\  —  rs)  is  equivalent  to  the  sheath  thickness  and  is  also  the  size  of  the  region  from  which  particles 
can  be  drawn.  Total  current  density  striking  the  debris  surface  is  given  as, 


(32) 

7(1  2?V(2s)) 

(33) 

/(V)  =  y„(i-2^.)) 

(34) 

mv  q 

7(F)  =  7o(1-A) 

(35) 

where  I  is  the  total  current  to  spherical  debris  equal  to  the  ambient  current  that  would  pass  through 
an  area  4kR2v  Jq  is  ambient  current  density  outside  the  sheath  and  equals  to  •  -A)  is  given  by 

Equation  (26).  The  current  derived  in  this  section  is  called  “thick-sheath,  orbit-limited”  current 
relation. 

In  the  derivation  of  Equation  (35),  velocity  of  the  body  is  not  taken  into  consideration  i.e., 
it  is  assumed  that  the  body  is  at  rest  compared  to  plasma  velocities.  At  geosynchronous 
altitudes,  plasma  energy  is  high.  The  average  kinetic  energy  of  particles  is  much  larger  than  the 
orbital  speed  and  the  body  is  essentially  at  rest  compared  to  plasma  particles  and  hence  Equation 
(35)  is  valid. 

At  low  Earth  orbit  regime,  body  gets  negatively  charged  i.e.,  ion  is  the  attracted  species.  However, 
at  low  Earth  orbit,  the  plasma  sheath  is  thin  such  that  thick  sheath  theories  do  not  apply  here.  Also, 
the  orbital  velocity  is  larger  than  the  average  kinetic  energy  of  ions.  Hence,  Equation  (35) 
becomes  invalid  for  ions  in  low  Earth  orbits. 
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For  low  Earth  orbits,  as  ions  are  relatively  at  rest,  they  are  assumed  to  be  ramming  into  debris 
surface  and  the  ion  current  is  approximately  given  as, 

I  =  emVscAi  (36) 

where  e  is  elementary  charge,  n,-  is  free  plasma  ion  density,  Vsc  is  inertial  debris  velocity  and  A,  is 
projected  ion  collection  area.  For  sphere  of  radius  Rs, 

A  =  (j)(4. *R])  (37) 

The  plasma  electron  and  ion  currents  to  a  sphere  of  radius  Rs  are  summarized  in  Table  1. 


Table  1.  Table  Summarizing  Plasma  Currents 


Orbit  Region 

Plasma  Current 

Plasma  Current 

(Species  is  Repulsed) 

(Species  is  Attracted) 

LEO 


(4^2) 


qnVsc(nR j) 


GEO 


[AnR]) 


(4  nR2s) 


It  is  to  be  noted  in  Table  1  that  a  debris  in  LEO  always  gets  negatively  charged  and  hence  the  ion 
automatically  qualifies  as  attracted  species. 

A  quick  and  approximate  method  to  validate  charging  level  in  low  Earth  orbits  is  sometimes  given 
by  the  so-called  Anderson’s  voltage  expression.  According  to  Philip  C.  Anderson  [6],  at  equilib¬ 
rium,  voltage  expression  is  given  as, 


Veq  = 


(38) 


where  K  is  Boltzmann  constant,  Te  is  plasma  temperature,  e  is  elementary  charge,  Ae  is  current 
collection  area  for  electrons,  Aj  is  current  collection  area  for  ions,  me  is  mass  of  electron  and  vsp 
is  speed  of  debris.  Given  below  is  a  quick  derivation  of  Anderson’s  formula. 


Anderson  approximated  plasma  ion  current  as, 

Ji  =  enjVsp  (39) 
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where  /,  is  ion  current  density,  e  is  elementary  charge,  nl  is  ion  density  and  vsp  is  debris  speed. 
Here,  it  is  assumed  that  ion  is  having  unit  elementary  charge  and  ions  are  relatively  stationary  w.r.t 
debris.  Assuming  the  presence  of  only  electron  and  ion  plasma  currents,  at  equilibrium. 


Electron  current  —  Ion  current 


(40) 


JeA 


e  — 


eVeq 

e  kr  j0Ae 


JA 

=  Mi 


1  2  kTe  eJeq 

-ene\  - e  k'r  Ae  =  enivspAi 

2  V  itme 


eVeq 

e  W 


Aj  „  /  Kme 


veq  = 


(assuming  ne  =  n,- ) 


(41) 

(42) 

(43) 

(44) 

(45) 


Equation  (38)  should  be  used  with  caution  for  it  is  valid  only  in  ionosphere,  and  that  too  approxi¬ 
mately.  But  Anderson’s  formula  acts  as  good  means  to  roughly  validate  results  for  low  Earth  orbit 
propagations,  specifically  when  debris  is  in  shadow.  Approximating  a  geosynchronous  orbit  with 
Anderson’s  formula  would  result  in  highly  inaccurate  results  because  Anderson  neglected  photo¬ 
electric  current  and  secondary  electron  current,  which  are  significant  in  GEO.  Table  2  argues  out 
the  reasons  for  exercising  restraint  in  use  of  Equation  (38)  in  full  detail. 
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Table  2.  Arguments  Against  Anderson’s  Voltage  Formulation 


Assumptions 

Arguments 

Ion  density  (n,)  =  electron  density  (ne) 

The  more  collect  assumption  should  be  neutrality  of 
charge,  so  n,  =  ne  assumption  should  be  valid  only 
when  charging  scenario  involves  electrons  and  a  sin¬ 
gle  monovalent  ion  species.  In  reality,  there  arc  more 
than  one  kind  of  ion  species,  and  they  may  be  multi¬ 
valent 

Ions  arc  essentially  stationary  w.r.t  debris. 
Mathematically, 

V particles / sc  =  Vparticels  ~  Esc  ~  Vsc 

Since  magnitude  of  ion  velocity  is  usually  of  the  order 
of  1.2  Km/s  or  less  in  LEO,  whereas  typical  spacecraft 
velocity  magnitude  in  low  Earth  orbit  can  be  7.8  Km/s, 
this  assumption  will  hold  quite  good  in  LEO.  This  as¬ 
sumption,  however,  will  not  hold  good  for  geostation¬ 
ary  orbit,  because  in  geostationary  orbit,  spacecraft 
velocity  magnitude  is  about  3  Km/s  and  ion  velocity 
magnitude  can  be  much  more  than  1.2  Km/s  as  there 
would  be  more  high  energy  ions  in  GEO  as  compared 
to  LEO 

Body  has  negative  potential 

This  assumption  again  holds  good  for  low  altitudes. 
At  high  altitude,  photoelectron  current  often  exceeds 
plasma  currents,  resulting  in  net  positively  charged 
body 

Effect  of  negative  potential  on  ion  collection  is 
neglected 

This  assumption  would  not  be  so  accurate  because  it 
can  be  seen  from  Whipple  [28]  that  body  potential 
indeed  effects  ion  current,  especially  the  accuracy  of 
this  assumption  would  be  poor  in  GEO 

Ignored  secondary  current,  backscattered  current  and 
photoelectric  current 

These  currents  cannot  be  neglected  as  insignificant, 
especially  the  photoelectric  current  in  GEO 

3.6  Secondary  Electron  Current  Model 


Expression  for  secondary  electron  current  can  be  found  out  using  yield  i.e.,  number  of  secondary 
electrons  generated  per  incident  projectile  particle.  Yield  expressions  can  be  multiplied  with  in¬ 
cident  plasma  ion  and  electron  currents  to  find  the  value  of  secondary  electron  current.  The  sec¬ 
ondary  electron  currents  due  to  hitting  ions  and  electrons  are  given  as, 


(<  A i  >)/,• 

(46) 

Ae  > )/e 

(47) 
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where  the  subscript  i  and  e  denote  hitting  particles  as  ions  and  electrons,  respectively.  /,  and 
Ie  represent  plasma  ion  and  electron  currents,  respectively  and  their  expressions  can  be  obtained 
from  Table  1 .  It  will  be  assumed  that  all  the  incident  electrons  and  ions  get  absorbed  and  contribute 
to  secondary  electron  emission  process.  The  <>  symbol  signifies  averaged  yield  value,  averaged 
over  all  energies,  and  is  given  as, 


si;mm(E)dE 

Si;EmdE 


(48) 


where  E  represents  energy,  f(E )  represents  the  Maxwell-Boltzmann  distribution  function  repre¬ 
senting  incoming  electrons  or  ions.  The  integral  limits  £)  and  Eu  in  Equation  (48)  represent  lower  and 
upper  energy  bounds  of  the  impacting  particles.  In  absence  of  precise  information,  Ej  and  Eu  are  taken  as 
0  and  °o,  respectively.  A (E)  represents  the  secondary  electron  yield  corresponding  to  the  particular 
energy  E  of  the  impacting  particle.  f(E)  is  given  as, 


f(E)=n 


m 

2nkT 


2  _  JL 
e  kT 


(49) 


It  may  be  argued  that  in  LEO,  ions  simply  ram  into  the  objects,  thereby  rendering  the  incoming 
ions  far  from  following  Maxwell-Boltzmann  distribution.  However,  the  energies  of  ions  and  elec¬ 
trons  in  LEO  are  too  small  to  incite  any  secondary  emission.  The  minimum  energy  requirement 
for  secondary  emission  can  be  computed  using  A (E). 


The  secondary  electron  yield  derived  is  taken  from  Stemglass  [29].  At  depth  x  of  the  substrate, 
low-energy  secondary  electrons  are  generated  in  two  major  ways  -  1)  produced  directly  by  inter¬ 
action  with  incoming  projectile  2)  produced  by  fast  8  rays.  8  rays  are  energetic  electrons  (induced 
by  energy  transfer  from  the  primary  projectile)  that  go  and  hit  other  electrons  and  atoms  to  form 
secondaries. 

Let  n^e  and  denote  the  number  of  secondary  electrons  produced  per  unit  distance  at  depth 
x  due  to  direct  interaction  and  fast  8  rays  respectively.  Let  Eq  be  the  mean  energy  loss  per  sec¬ 
ondary  formed.  Let  be  the  mean  energy  loss  per  unit  distance  going  directly  into  secondary 

formation.  Therefore, 


nsl\vi,x)  =  ^r{^-)S  (50) 

Lq  ax 

n£\vi,x)  =/(v/,x)i(^)i?  (51) 

where  is  the  energy  loss  per  unit  length  going  into  the  formation  of  8  rays,  and  f(vi-x)  is  a 

factor  that  represents  the  fraction  of  available  for  formation  of  secondaries  in  higher  order 

collisions  at  depth  x.  According  to  Stemglass, 


,dEi_\(\)  _  ,dEi\{ 2)  _  1  ,dEj 
'  dx  'av  ~  '  dx  )av  ~  2 {  dx 


(52) 
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(53) 


where  (^)av  is  the  total  energy  loss  per  unit  length. 


nse(v:x)  =  (i)(^)(^)av[!  +f(vi,x)] 


According  to  Bohr, 


,dEj 


)av 


2nN£4zf 

p 

*^eq 


(54) 


where, 

r  1  2  fm0\T7 

Eeq= -niQVj  ={—)Ei  (55) 

N  is  the  number  of  atoms  per  unit  volume,  e  is  electronic  charge,  z.j  is  charge  of  incident  particle,  M 
is  mass  of  incident  particle,  Ej  is  energy  of  incident  particle,  vl  is  velocity  of  incident  particle,  z„j 
is  number  of  electrons  in  n,l  shell  and  /„  /  is  binding  energy  of  these  electrons.  The  summation  is  to 
be  taken  over  all  shells  for  which  logarithm  remains  positive.  If  the  primary  velocity  is  sufficiently 
high,  then  Equation  (54)  takes  the  following  form, 


,dEj 
'  dx 


)av 


iKNe^Elnd^l)] 

t-'eq  * 


(56) 


where  /  is  mean  excitation  potential  for  the  atom,  Z  is  atomic  number  of  substrate  atom  and 
Zln(T)  =  £n jZn,iln(In,l)-  E()  is  approximated  as  25eV  for  solids.  The  idea  of  minimum  energy 

required  for  secondary  electron  emission  stems  from  Equation  (56).  In  order  to  ensure  that  ‘ 
)avis  positive,  input  inside  the  natural  logarithm  function  in  Equation  (56)  must  be  greater  than  1 
words,  Eeq  must  be  greater  than 


A  secondary  electron  formed  at  depth  x  below  the  surface  will  have  certain  probability  P(x)  of 
reaching  the  substrate  surface  and  escaping.  P(x)  can  be  given  as, 

P(x)  =  T Ae~t  (57) 


where  r  is  surface  transmission  coefficient,  representing  the  probability  that  an  electron  arriving  at 
the  surface  from  the  interior  will  be  able  to  escape.  A  is  a  constant  determined  by  the  distribution  of 
initial  velocities  of  secondaries  and  by  the  ratio  =  nc,  where  Xsa  is  mean  free  path  for  absorption, 
Xsc  is  mean  free  path  for  inelastic  collisions.  For  a  symmetrical  distribution  of  initial  directions 
about  a  plane  parallel  to  the  surface,  A  has  an  approximate  value  of  .6.  Ls  appearing  in  Equation 
(57)  is  given  as, 

Ls  =  C-lsa^scP  (58) 

Ls={l-nc)hsc  (59) 

The  mean  free  path  for  inelastic  collisions,  Xsc  is  given  as, 

Kc  =  (Nasc)~l  (60) 

where  osc  is  the  cross-section  for  secondary  electron  scattering,  and  is  given  as, 
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Osc  =  CCGg  (61) 

where  og  is  geometrical  area  of  the  outermost  filled  electron  shells  as  determined  by  covalent  radii 
and  a  is  some  constant.  <Jg  is  approximately  given  as, 

og  =  1.6Z*  x  10' 16  (62) 

where  Z  represents  atomic  number  of  the  target.  Using  Equations  (60),  (61)  in  Equation  (59), 


Ls  =  (^nc)^(Naag)  1  =  (a'Nog)  1  (63) 

where  a1  =  (|  )2a.  If  the  assumption  that  the  solid  may  be  treated  as  a  collection  of  free  atoms 
is  valid,  d  should  have  closely  the  same  values  for  all  metals,  a' is  assumed  to  be  .23  [29].  Surface 
transmission  coefficient  T  appearing  in  Equation  (57),  in  case  of  isotropic  distribution  of  electron 
velocities  approaching  a  uniform  surface  potential  barrier,  is  given  as, 


T  =  1  — 


2 

Es  +  </>D  / 


(64) 


where  Es  is  energy  of  the  secondary  as  measured  relative  to  zero  vacuum  level  outside  the  metal. 
Typical  values  of  Es  lies  in  the  range  of  6  eV  to  8  eV.  Also,  for  monovalent  metals,  typical  values 
of  (p/j  lies  in  the  range  of  .1  eV  and  .5  eV.  Thus,  it  can  be  seen  that  r  lies  in  the  range  of  .8  and  .9. 
This  makes  tA  «  .5,  and  hence, 

P{x)  ~  ,5e~E  (65) 

Yield,  A  is  defined  as  as  the  number  of  secondary  electrons  escaping  from  the  surface  per  incident 
ion.  The  yield  from  a  thin  layer  of  width  dx  located  at  depth  x  is  given  as, 


dA  =  nse(vi,x)P(x)dx  (66) 

r°°  1  1  AF  x 

A=  /  xiri^Ull+fMjTAe-rsdx  (67) 

Jo  2  Eq  dx 

When  the  atomic  number  of  substrate  is  small  (Z  <  30), 

f(Vi,x)  =  l-e  Os  A)  (68) 

where  L§  (v;)  is  effective  penetration  distance  of  8  rays.  Substituting  Equation  (68)  into  Equation 

(67),  and  recalling  that  (~d~xL  )av  maY  be  taken  as  constant  over  a  surface  region  of  width  large 
compared  to  Ls, 

1  1  dE 

A  =  0^(-7i)AvTAL,(l  +F(Vi))  (69) 

z  Lq  ax 

^(v/Ha+^r1  (70) 

For  heavier  elements  (Z  >  30),  more  accurate  expression  for  /(v,-,x)  results  in, 

A=  T7T<^>^tAL*(1+^(v')) 
z  Lq  ax 
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(72) 


f(>'i)  =  ((i+^r1+^)(i+^)- 

As  Eg  Lg 


It  has  been  found  that  [29], 


(vi)  ^  Eeq 

L,  ~  100 


(73) 


where  Eeq  is  expressed  in  electronvolts.  Experimental  values  for  the  ratio  range  typically  from 
.05  for  aluminium  to  .25  for  gold  [29].  A  linear  variation  of  this  quantity  with  atomic  number  will 
be  assumed,  i.e., 

Lx/ 

-r-  =  -05  + 

L8 


.25 -.05 
79-13 


(Z—  13) 


(74) 


It  is  to  be  noted  carefully  that  the  formulations  are  made  assuming  centimeter-gram-second  (CGS) 
system  of  units.  Appropriate  unit  conversions  must  be  applied  so  as  not  to  run  into  unnecessary 
errors. 


3.7  Backscattering  Current  Model 

Backscattering  current  can  be  obtained  with  the  help  of  backscattering  yield  i.e.  number  of  backscat- 
tered  particles  per  incident  particle.  The  backscattering  of  electrons  is  more  important  than  the 
backscattering  of  ions.  Backscattering  yield  is  usually  much  smaller  compared  to  secondary  yield. 
Backscattering  electron  current  is  given  as, 


1 bsc  —  t'Ie  (75) 

where  r  represents  backscattering  yield  and  Ie  represents  plasma  electron  current,  which  can  be  ob¬ 
tained  from  Table  1 .  One  might  be  wondering  why  averaged  value  is  not  considered  for  backscat¬ 
tering  yield.  The  fact  is,  as  will  be  seen  in  next  section,  according  to  the  model  chosen,  backscat¬ 
tering  yield  comes  out  to  be  constant  for  a  particular  material  i.e.,  <  r  >=  r. 

There  are  two  models  which  Archard  [30]  used  for  explaining  backscattering.  Neither  of  these 
two  models  single-handedly  explains  the  backscattering,  but  a  combination  of  the  two  explains  the 
backscattering  over  entire  atomic  number  range. 

In  the  diffusion  model,  it  is  assumed  that  electrons  travel  straight  into  the  target  up  to  a  certain 
specified  distance,  after  which  they  diffuse  evenly  in  all  directions.  This  model  assumes  that  an 
electron  exhibits  completely  random  motion  after  multiple  collisions.  It  ignores  the  possibility  of 
an  electron  undergoing  large  single  elastic  reflections  between  the  surface  and  the  depth  of  com¬ 
plete  diffusion. 

In  the  elastic  collision  model  [31],  it  is  assumed  that  electrons  travel  straight  into  the  target  suffer¬ 
ing  retardation  according  to  Thomson- Whiddington  law.  In  addition,  the  electrons  undergo  elastic 
collisions  in  accordance  with  Rutherford’s  law  of  scattering.  This  model  acknowledges  the  pres¬ 
ence  of  electrons  that  are  elastically  scattered  through  large  angles,  but  ignores  the  diffusion  effect 
of  multiple  collisions.  All  electrons  that  are  not  turned  through  more  than  a  right  angle  are  assumed 
not  to  have  turned  at  all. 
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“Depth  of  complete  diffusion”  is  the  depth  at  which  the  average  cosine  between  actual  direction  of 
motion  and  direction  of  primary  beam  becomes  2  or,  the  angle  becomes  68.4°.  Corresponding  to 
this,  a  critical  energy  Eci  exists  such  that, 


fE o  dE  _  1 

4  A|f  |  2 


(76) 


1 

I 


KNe4Z2 
2  E2 


ln 


2 fWjf) 

HZ? 


1 

2 


(77) 


where  N  is  number  of  atoms  per  cubic  centimeter,  e  is  electronic  charge,  m  is  electronic  mass, 
h  —  A  is  reduced  Plank’s  constant,  Z  is  atomic  number  of  target,  an  is  Bohr  hydrogen  radius 
=  5.2917721092  x  10~n  m,  Eq  is  initial  energy. 


From  Bethe’s  stopping  power  theory, 


dx  E 


(78) 


where  J  =  11. 5Z  (eV)  is  approximate  mean  excitation  energy.  Substitution  of  Equations  (77), 
(78)  into  Equation  (76)  results  in, 


fEo  nNe4Z2  /  2aHm(  ^ )  2  \  EdE  _1 
JEd  2 E2  "V  Kzh  )  2nNe4ZlnCf)  ~  2 


Applying  change  of  variables  V  —  f  from  Energy  E  to  voltage  V, 

fvo  Z  (eV)edV  /2aHm(^-)^\  eVedV  _1 

Jvd  2(eV)22 feg)  A  nzl  )  2 KNe4Zln(^)  ~  2 


fvo  Z  ln(. 5418E2Z31)  1 

4V  /«(.  1739VZ  ')  “20 


(81) 


The  ratio  of  logarithms  in  Equation  (81),  according  to  Archard  [30],  remains  constant  between  5 
KV  and  50  KV  and  is  approximately  taken  by  Archard  as  .7.  However,  this  does  not  always  hold 
true.  As  a  better  approximation  to  this,  it  can  be  assumed  that  the  logarithmic  ratio  is  .7  for  energy 
between  5kV  and  100  kV,  and  the  ratio  is  .6  for  higher  energies. 


Proceeding  forward  with  the  assumption  of  ratio  as  .7  (.6  assumption  follows  the  same  procedure), 
Equation  (81)  results  in, 

( V{)  \  20 

Z//7(— )  -  —  (82) 

On  substituting  V  =  E  in  Bethe’s  stopping  power  formula, 


dE 

dx 


27lNeAZ 


E 


(83) 
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(84) 


dV 

dx 

Integration  of  Equation  (84)  leads  to, 


2 nNeiZ1  ,  2V  , 

— —,n(Tnz) 


X  = 


7.68  x  1012 
NZ 


(V02F(F0)-F2F(V)) 


(85) 


where  x  is  in  centimeters,  N  is  number  of  substrate  atoms  per  cubic  centimeter,  Vo  is  initial  voltage 
in  volts  and  the  function  F  is  defined  as  follows, 


Now, 


i  i  ?' 

F(V)  =  -(1  +  -  +  — + . ) 

y  y  y~ 

(86) 

,i74yN 

y-2ln{  z  ) 

(87) 

N=pN* 

A 

(88) 

where  p  is  target  density  in  gm/cc,  Na  is  Avogadro’s  number,  A  is  atomic  weight  in  gram/mole. 
Substituting  Equation  (88)  in  Equation  (85)  results  in, 


px 


7.68  x  1012  x  A 
NAZ 


(V02F(F„)_F2F(V)) 


(89) 


Assuming  electrons  come  to  rest  when  F(V)=0, 


Pxr 


7.68  x  1012V02  xA 
AAZ 


F(V0) 


(90) 


where,  xr  denotes  full  range.  The  ratio  of  depth  of  complete  diffusion  to  full  range  can  be  written 


as, 


As  per  Archard, 


Xd  =  j 

xR 


xd  V$F{Vo)-V][F{Vd) 
xR  V02E(y0) 

^-.1  (vA2m) 

XR  \Vo)  F(Vo) 


Vd 

Vo 


2  ln(FpL 


2  ln(FpL 


4  + 


(1  +  1  +  5  +  - 


2! 


2  ln{ 


.174  V), 


VO  V5 


-(1  +  f  +  %  +  ...) 

)  y  yd  yd 


+  i  +  + 


2  ln( 


Using  the  approximation  given  by  Equation  (94)  in  Equation  (93), 


^~1  _(WA 


xr 


Vo 


(91) 

(92) 

(93) 

(94) 

(95) 
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Xd 

XR 


(96) 


—  l—e  7Z  [since  Zln{^~)  = 


xR  V  V  1Z 


[first  order  approximation ) 

Xd_  _  40 

7Z 


(97) 

(98) 


With  the  assumption  that  electrons  move  equally  in  all  directions  from  the  depth  of  complete 
diffusion  in  such  a  manner  that  the  overall  path  length  equals  full  range,  it  can  be  easily  seen  from 
Figure  3  that  only  electrons  within  solid  angle  20  will  escape. 


Figure  3.  Backscatter  Cone  [30] 

The  backscattering  coefficient/yield  R  is  given  as, 


Solid  angle  formed  by  semiangle  0 

R  = 

total  solid  angle 

(99) 

27r(l—  COS0) 

471 

(100) 

1  Lx) 

«=  2 

(101) 

7Z-80  xd  40 

14Z-80  V  xR  1Z ’ 

(102) 

From  Equation  (102),  it  can  be  seen  that  R  is  zero  for  Z  <  12  (as  R  cannot  be  negative).  In  fact,  for 
low  Z  (Z  <  1 1),  experimental  value  does  not  match  with  theoretical  value  zero.  This  is  because,  as 

Z  decreases,  the  ratio  x~dR  increases  i.e.,  electrons  travel  a  large  part  of  their  range  before 
diffusion  sets  in.  Hence,  “there  is  significant  time  for  other  neglected  phenomenon  of  large 
single  elastic  collisions  to  take  place. 
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Backscattering  Yield  from  Elastic  Collision  Model  According  to  Thomson-Whiddington  law, 
if  v  is  the  electron  velocity  and  x  is  the  distance  traveled  in  the  target,  then, 

v4  =  Vq  —  cpx  —  cp(R  —  x)  (103) 


where  v  is  electron  velocity  in  cm/s,  vo  is  initial  electron  velocity  in  cm/s,  p  is  target  density  in 
gm/cc,  c  is  a  constant  and  its  value  is  5.05  x  1042cm6sec4g~1,  R  is  range  of  electron  in  target 

v4 

material  in  cm  and  given  as  R  —  Rutherford  [31]  predicts  that  if  no  electrons  per  S  square 
centimeter  are  incident  on  a  bare  nucleus  of  charge  Ze,  the  scattered  intensity  per  unit  solid  angle 
is  given  as, 


(  n0Z2e4\ 


m  WifivJLn4^) 


(104) 


where  (j)  is  angle  of  deviation  of  scattered  electrons  as  shown  in  Figure  4,  e  is  electronic  charge,  no 
is  electron  density  per  S  sq.  cm,  m  is  electronic  mass.  If  a  rectangular  parallelepiped  element 
having  surface  area  S  and  thickness  dx  is  considered,  then  the  number  of  atoms  present  in  the 
volume  element  is  given  as, 


dN  =  Na  ■  moles 

NApSdx 

dN  — - 

A 


(105) 

(106) 


where  NA  is  Avogadro’s  number  and  A  is  gram  atomic  weight  of  the  substrate  material.  Incremental 
number  of  electrons  at  depth  y  deflected  through  an  angle  of  (j)  =  K  —  6  into  the  incremental  solid 
angle  dQ.  is  given  as, 

dn(y,e)=P(^)-dN-dQ  (107) 


Figure  4.  Backscatter  Cone,  scattered  electrons  [31] 


Substituting  the  values  of  P{<j>),  dN  and  Q.  =  27t(1  —  cos0  )  or  d£l  =  2k  sin0  dO  in  Equation  (107), 

noZ2e4  1  NApSdxr 


dn(y-Q)  = 


ASm2cp  (R  —  x)  sin4  ^  A 


-2KsinOdO 


(108) 
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On  change  of  variable  using  y  =  f ,  Equation  (108)  results  in, 

Z2eAN/\n{)(y)R(dy)2nsm  9d9 
^  ’  =  4m2cA(  1  —  y)  cos4  ( |) 

Number  of  electrons  incident  on  a  plane  at  depth  y  in  the  target  is  given  as, 

ry  /■! 


(109) 


no(;y)  =  «o 


Jo  jo 


dn(y,9) 


(HO) 


It  is  assumed  here  that  electrons  having  deflection  angle  of  less  than  90°  have  not  deflected  at  all, 
and  therefore,  reflection  coefficient  r  obtained  is  smaller  than  experimental  value.  Using  Equation 
(109)  in  Equation  (110)  results  in, 

/  \  fyno(y ) 

no(y)  =  no  —  a  /  - - dy 

Jo  1-y 


(111) 


where, 


KZ2eANA 
m2cA 


a  = 


Equation  (111)  has  the  following  solution, 

«o(y)  =  «o(i  - y)a 

Incremental  backscattering  coefficient/yield  is  given  as, 

dn(y,  0) 


dr(y,  9)  = 


no 

\a—  1 . 


,  /  an  \a-lj  sin 9d9 

drM  =  -2(l-y)  dy _ 


(112) 


(113) 


(114) 


(115) 


0is  varied  from  0  to  0o  where  x  +  R  ory(l  +  sec0o)  =  1,  and  y  is  varied  from  0  to  .5  in 

integration  of  Equation  (115).  The  integration  results  in  following  expression  for  backscattering 
yield. 


Backscattering  Yield, r  = 


a— l  +  .5fi 
a  T  1 


(116) 


The  value  of  constant  ‘a’  as  given  by  comes  out  to  be  .012.  However,  Everhart  [31]  later 

changed  it  to  .045  to  match  the  experimental  results.  Archard  has  given  a  very  nice  explanation 
for  this  whimsical  change  of  ‘a’  by  Everhart.  According  to  Archard,  half  the  electrons  deflected 
between  0°  and  90°  will  be  lost  in  layers.  Let  us  suppose  an  electron  suffer  45°  deflection,  it 
should  suffer  another  45°  deflection  in  opposite  sense  to  retain  original  direction.  It  may  also  suffer 
deflection  of  45°  in  same  sense  and  get  lost.  So,  Equation  (1 10)  would  actually  have  the  following 
form, 

ry  r?  l  /■>’  /■? 

no{y)  —  no—  /  dn(y,  9)  —  -  /  dn{y,9)d9  (117) 

Jo  Jo  2  Jo  J f 

If  the  above  formulation  is  used,  ‘a’  comes  out  to  be  .041,  which  is  very  close  to  the  empirical  ‘a’ 
of  .045  used  by  Everhart. 
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At  low  atomic  number  (Z),  the  ratio  of  depth  of  complete  diffusion  xd  to  full  range  xr  is  large, 
such  that  there  is  a  great  chance  of  large  elastic  collisions  before  diffusion  sets  in,  hence  elastic 
collision  model  is  used.  For  high  Z,  ^  is  small,  so  diffusion  sets  in  almost  immediately  long  before 
any  elastic  collision  may  take  place,  so  diffusion  model  is  applicable.  In  the  intermediate  region, 
both  models  are  applicable  and  a  weighted  average  of  both  is  used.  The  following  piece-wise 
model  can  be  assumed  for  a  substrate  with  atomic  number  Z: 

For  Z  <  11, 


For  Z  >  60, 


For  1 1  <  Z  <  60 


r 


r  — 


a  —  1  +  .5a 

ci  1 

7Z-80 

14Z-80 


r_x±f  a~  1  +  -5a\  xR-xd  (  7Z-80 
Xr  l  a+  1  )  Xr  \14Z  —  80 


(118) 

(119) 

(120) 


3.8  Photoelectron  Current  Model 


Photoelectron  current  is  one  of  the  significant  currents  in  high  altitude  orbits  like  GEO.  Photoelec¬ 
tron  current  density  (divided  by  elementary  charge)  is  given  as  [7], 


f‘00 

Jph  =  Jo  fsYdco  (121) 

where  fs  is  number  of  photons  per  unit  area  per  unit  time  per  unit  photon  energy,  and  is  a  function 
of  photon  frequency,  co.  Y  is  the  yield  i.e.  number  of  photoelectrons  per  incident  photon  and  is  a 
function  of  co  and  angle  of  incidence,  0.  The  yield  is  given  as 

y[Q), /?(©)]  =y*[Q),/?(o),0)][l-/?(Q),0)]  (122) 

where  R  is  reflectance  and  Y*  is  yield  per  absorbed  photon,  which  can  be  approximated  as  [7,  32], 


y*[G>,/?(a>,0)]  ~ 


y*[G>,/?(G>,o)] 

cosO 


(123) 


Laboratory  measurements  [7,  33,  34]  show, 


1 -/?(©,  0)  «  [1 -/?(©, O)]cos(0)  (124) 

Thus,  from  the  product  of  Y *  and  (1  —  R),  there  is  no  0  dependence  of  Y.  The  only  incident  angle 
dependence  of  the  photoelectric  current  stems  from  the  factor  cos  0,  which  results  from  effective 
surface  area  on  which  light  is  incident.  Photoelectrons  can  be  assumed  to  follow  Maxwellian 
distribution  for  positive  body  potentials  with  characteristic  temperature  ranging  leV-2eV  [7].  The 
characteristic  temperature  Tp]z  will  be  taken  as  1.5  eV  in  our  simulations.  This  leads  to: 


Iph  —  qeJph(0)Acos9  —  (]eJph(Q)Aj_  if  0  ^  0  (125) 
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(126) 


Iph  qeJph(0)Aj_e  ph  if  (j)  >  0 

where  the  neutral  current  density  Jph(0)  is  given  as, 


Jph( 0)  =  J  fs(co)Y*(co)[l-R((o)]d(o  (127) 

where  qe  represents  elementary  charge,  A_  represents  projected  area  of  the  body  and  <p  represents 
the  body  potential. 

In  the  following,  methods  for  obtaining  quantities  appearing  in  Equation  (127)  are  stated.  Number 
of  photons  per  unit  area  per  unit  time  per  unit  photon  energy  ( fs )  will  be  obtained  from  Air  Mass 
Zero  (AM  0)  solar  spectral  irradiance.  Yield  per  absorbed  photon  (Y  *)  will  be  obtained  using 
Fowlers  method  for  near- threshold  frequencies  and  assumed  constant  for  higher  frequencies. 
Fowlers  yield  depends  upon  absolute  body  temperature,  which  will  be  obtained  using  thermal 
equilibrium  equa-tion.  Reflectance  (R)  will  be  obtained  using  Fresnel’s  equation.  Fresnel’s 
equation  depends  upon  refractive  index,  which  will  be  computed  using  Brendel  and  Bormann 
model. 

AM  0  or  ASTM  E-490  spectrum  provides  solar  spectral  irradiance  (  2  -i )  data  and  is  based 
on  data  collected  from  several  satellites,  space  shuttles,  high  altitude  aircrafts,  solar  telescopes  etc. 
The  correctness  of  the  spectrum  is  validated  by  scientific  communities  using  the  spectral  irradiance 
integration,  which  should  be  equal  to  the  accepted  value  of  solar  constant  i.e.  1366.1  J^2-  AM 
0  spectrum  is  shown  in  Figure  5. 


(a)  AM  0  Solar  Spectrum  Irradiance  (b)  AM  0  Solar  Photon  Flux 

Figure  5.  AM  0  Spectrum 


If  A,  and  A  f  are  the  initial  and  final  wavelengths  of  a  certain  band,  then  the  power  received  per  unit 
area  can  be  written  as, 


(128) 
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where  S-A  is  spectral  irradiance  (power  per  unit  area  per  unit  wavelength).  Then,  photon  rate  per 
unit  area  (number  of  incident  photons  per  unit  area  per  unit  time)  in  that  A  range  can  be  written  as, 


Ph 


rate 


(129) 


where  E  is  energy  of  a  photon  having  wavelength  A .  The  photon  rate  per  unit  area  can  also  be 
written  in  the  following  form, 

rEf 


Phrate  =  -  /  fsdE 


(130) 


where  fs  represents  solar  photon  flux  (number  of  photons  per  unit  area  per  unit  time  per  unit  photon 
energy).  The  negative  sign  is  because  of  the  fact  that  Ay  >  A;  =>  Ef  <Ej  and  vice  versa.  Bringing 
in  a  change  of  variables, 


(131) 


Comparison  of  Equation  (129)  and  Equation  (131) 

results  in,  yr 

Sx  =  -fsE—  (132) 


Differentiating  the  energy  equation  E  = 
(132)  results  in, 


Sx=fs 


he 
X  ■ 


dE 

dX 


_  hc- 


AS 


(hc): 

V 


or 


fs  = 


obtained,  and  substituting  it  in  Equation 
(he) (133) 


Application  of  Equation  (133)  to  solar  spectrum  irradiance  data  presented  in  Figure  5a  results  in 
the  solar  photon  flux  data  plot  given  by  Figure  5b. 


In  the  field  of  quantum  statistics,  Fermi-Dirac  distribution  states  that  for  a  system  of  identical 
fermions,  the  average  number  of  fermions  in  a  single-particle  state  i  is  given  as, 


1 

'l  ~  e(Ei~H)/{kT)  _|_  1 


(134) 


where  e,  is  the  energy  of  the  single -particle  state  i,  /i  is  total  chemical  potential,  k  is  Boltzmann 
constant  and  T  is  absolute  temperature.  The  term  fermion  is  used  for  any  particle  following  Fermi- 
Dirac  distribution.  If  it  is  assumed  that  the  electrons  in  a  material  follow  the  Femi-Dirac  distribu¬ 
tion,  then  the  number  of  electrons  per  unit  volume  with  velocity  components  in  the  ranges  (u,u+du), 
(v,v+dv),  (w,w+dw)  (u  being  normal  to  the  surface  and  u,v,w  forming  an  orthogonal  system)  can 
be  written  as, 


«(m,v,  w)dudvdw  —  2  (  — 


m 


dudvdw 


(135) 


e(jm(u2+v2+w2)—e*)/(kT)  _|_  i 

Here,  m  is  mass  of  electron,  h  is  Plank’s  constant  and  e*  is  fermi  energy,  k  is  Boltzmann  constant 
and  T  is  absolute  temperature.  The  particle  volumetric  density  having  surface-normal  velocities  in 
the  range  (u,u+du)  is  given  as, 


n(u)du 


if—)3 du  r r  ■  ,pipie — 

V  h  )  Jo  Jo  e{\m{u‘-+p2)-£*)/{kT)  +  j 


(136) 
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n(u)du=^tI^(t^]  log{l+e{£*--2'nu2)/^)du  (137) 

m  \h  J 

A  new  variable  Xo  is  introduced  such  that  Xo  —  X  +  £*,  where  %  is  the  work  function.  Fowler’s 
derivation  [35]  is  based  on  the  hypothesis  that  “the  photoelectron  yield  is  proportional  to  the  num¬ 
ber  of  electrons  per  unit  volume  of  the  metal  whose  kinetic  energy  normal  to  the  surface  augmented 
by  photon  energy  hv  is  sufficient  to  overcome  the  potential  step  at  the  surface”.  Number  of  such 
available  electrons  Nb  is  given  as, 


Nb  = 


n(u)du 


J  ^mu2=Xo—hv 

With  introduction  of  a  new  variable,  y  —  (\inp2  —  (xo  —  /?  v ) ) / (kT) 


(138) 


Nb  = 


2nkT  flkT\2  fm 


m  \  m 


r°  iog{\  +e-y+(hv-x)/(kT)) 

'°  c y  +  (Xo-hv)/(kT))i 


dy 


(139) 


When  the  frequency  is  near  threshold  frequency,  y  is  small  compared  to  the  term  (xo  ~  hv)/(kT ), 
so  Fowler  neglects  the  denominator  y  in  Equation  (139).  This  is  the  single  most  important 
assumption  of  Fowler’s  theory  and  it  leads  to: 


Nb  = 


2y/2nm2  k2T 2 


logil+e-y+^-xV^dy 


/?3  (Xo~hv)iJo 

The  logarithmic  term  may  be  expanded  as  a  series  and  integrated  term-wise  to  obtain, 


e2^  e^ 


2\p2% m2  k2T2  .  ..  ^ 

B=  a3  to-A vfe  -vr+^T~'"] 


where  p  =  (hv  —  x)/(kT).  When  p  >  0, 


r°°  rt1  i*oo 

/  log(\  +  e~y+^)dy  —  /  log(l + e~y+^)dy +  /  log(l  +e~y  )dy 

Jo  Jo  Jo 

r°°  jr-  1  rU  ,, 

J Q  log(l  +  e~y+^)dy= -  +  -p2  + log(l+e~y  )dy” 

The  logarithmic  term  is  expanded,  which  results  in, 


Jo 


log(  l  +  e 


7T-  1  p  -U  p 

-y^)dy=\  +  \p2-[e-^- 
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+ 


32 


Therefore, 


Nb  = 


2\f2n. 


3 

m2 


2rr2 


kzT 


K2  1  ?  _u  e  2/J  e 
h3  (Xo-hv)2l^  +  ^~-[e  + 


(140) 


(141) 


(142) 

(143) 

(144) 


(145) 


It  is  noteworthy  that  even  when  p  <  0,  photoemission  occurs,  because  Fowler  tries  to  evaluate  a 
situation  where  kinetic  energy  of  motion  perpendicular  to  surface  plus  incident  light  energy  hv 
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exceeds  the  work  function.  Since  Fowler  has  assumed  Np  being  proportional  to  the  yield,  it  can  be 
written. 


Y* 


RP - — 

(Xo-hv)i  \  kT  J 


(146) 


where  Kp  is  a  constant  of  proportionality.  Fowler  also  assumed  that  for  frequencies  near  threshold, 
(Xo  ~  hv)  2  is  constant.  According  to  Fowler,  at  hv  «  %,  (Xo  —  hv)  ~  £*,  so  if  v  changes  by  15%, 
(with  Vo  of  the  order  of  4  eV  and  Xo  ~  hv  of  the  order  of  10  volts),  (xo  ~  hv)  changes  by  6%  and 
(Xo~hv)2  changes  by  only  3%.  This  implies  approximately  a  constant  (Xo  —  hv)i.  Our  upper 
and  lower  limit  on  variation  of  photon  frequency  V  from  threshold  frequency  Vo  would  correspond 
to  3%  variation  of  (Xo  —  hv)  2  from  its  nominal  value  (Xo  —  x) 1  ■  Thus, 


r=«;r2/(ji) 


(147) 


where  K'p  is  new  constant  of  proportionality.  The  function  f  as  previously  given  is  defined  as, 

f(n)  =  eP  with  jd  <0  (148) 

and, 

7r2  1  s>  2/i  p  3  fi 

/(m)  =  -g-  +  2^ (e_/i  “  ^2“  +  ^2“  _  "  •)  with  M>0  (149) 

The  value  of  constant  K’p  has  to  be  found  from  experimental  values  of  the  photoelectric  yield.  For 
frequencies  that  have  energy  values  lower  than  the  proposed  3%  range,  it  will  be  assumed  that  there 
is  no  photoelectron  emission.  For  all  frequencies  that  have  energy  values  higher  than  this  range, 
the  yield  will  be  taken  as  constant.  This  assumption  is  made  using  the  basic  fact  of  photoelectric 
effect,  “for  frequencies  greater  than  threshold  frequency,  there  is  no  effect  of  frequency  in  photo¬ 
electric  current.  But  it  increases  linearly  with  increase  in  light  intensity.  Increasing  frequency 
only  leads  to  increase  in  kinetic  energy  of  photoelectrons” .  As  can  be  seen  in  Figure  5b,  above 
threshold  frequency  (4.16  eV  for  aluminum),  intensity  is  approximately  constant,  thereby 
validating  the  cor-rectness  of  assumption  of  constant  yield  above  the  proposed  3%  range. 


Equation  (147)  requires  absolute  body  temperature,  which  can  be  obtained  from  thermal  equilibrium. 

Heat  absorbed  =  (1  —a)7lR2So  (150) 

where  a  is  albedo,  R  is  radius  of  the  body  and  So  is  solar  constant 

Heat  radiated  from  the  body  —  (4nR2)aT4  (151) 

where  a  is  Stefan’s  constant.  Under  thermal  equilibrium  condition, 

(1  —  o)7iR2Sq  =  (4-tiR2)oT4  (152) 

(1  -a)S0\* 


T  = 


4(7 


(153) 
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Albedo  is  defined  as  fraction  of  incident  light  that  is  reflected  off  of  the  body  surface.  Reflectance 
is  defined  similarly  but  for  a  single  incidence  angle.  Thus,  albedo  is  directional  integration  of  re¬ 
flectance;  albedo  can  be  obtained  as  sum  of  diffuse  and  specular  reflection  components. 

An  important  property  closely  related  to  reflectance  is  reflectivity.  Reflectivity  is  property  of  a 
material  whereas  reflectance  is  property  of  a  particular  sample  of  that  material.  Only  in  case  of 
very  thin  bodies,  reflectance  and  reflectivity  are  different.  If  the  body  is  too  thin  e.g.,  of  the  order  of 
few  atoms  layer  thin,  inner  atomic  layers  start  reflecting  the  light  as  well  as  the  light  starts  transmit¬ 
ting  through  the  body  easily,  resulting  in  different  values  of  reflectance  and  reflectivity.  For  thick 
materials,  both  reflectance  and  reflectivity  are  same.  Reflectivity  is  limiting  value  of  reflectance 
when  your  sample  is  thick.  Reflectivity  is  given  as  square  of  magnitude  of  Fresnel’s  reflection 
coefficient.  It  can  be  safely  assumed  that  debris  is  thick  and  reflectance  can  be  substituted  with 
reflectivity  and  hence  Fresnel’s  equation  can  be  used. 

The  expression  for  Fresnel’s  equation  depends  upon  the  optical  properties  of  the  body  and  the 
medium  in  which  incident  light  is  present.  Although  current  research  deals  with  air  (dielectric) 

-  conductor  interface,  we  will  first  look  at  the  derivation  of  Fresnel’s  equation  for  dielectric  - 
dielectric  interface  and  see  how  it  relates  to  Fresnel’s  equation  for  dielectric  -  conductor  interface. 

Any  electric  field  can  be  separated  into  two  components,  one  perpendicular  to  the  plane  of  Figure 
(6)  and  also  perpendicular  to  the  direction  of  light  travel,  and  the  other,  parallel  to  the  plane  of 
Figure  (6)  and  also  perpendicular  to  the  direction  of  light  travel  [36].  Figure  (6a)  corresponds  to 
E  (E±  lies  out  of  the  plane  of  paper  and  corresponding  B  lies  in  the  plane  of  paper).  Similarly, 
Fig.  (6b)  corresponds  to  £ji  (£j|  lies  in  the  plane  of  paper  and  the  corresponding  B±  lies  out  of  the 
plane  of  paper).  The  small  lines  perpendicular  to  light  travel  directions  in  Figure  (6a)  correspond 
to  B-fields.  Similarly,  the  small  lines  perpendicular  to  light  travel  directions  in  Figure  (6b) 
correspond  to  E  fields.  All  the  in-plane  field  components  are  divided  into  their  horizontal  and 
vertical  components  as  can  be  seen  in  Figures  (6a),  (6b). 


(a)  E  j_,  B  ||  to  Plane  of  Paper  (b)  B  L,  E  to  Plane  of  Paper 
Figure  6.  Electric  Field  Break  Down  into  Parallel  and  Perpendicular  Directions 


Derivation  corresponding  to  Figure  (6a)  will  be  first  looked  at.  Since  the  E-field  is  parallel  to 
the  boundary  (boundary  is  out  of  the  plane  of  paper),  the  E-field  will  just  transmit  across  the 

boundary 
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(consequence  of  Maxwell’s  third  equation)  i.e., 


Ej  +  E,  —  Et 


(154) 


where  the  subscripts  i,  r  and  t  stand  for  incident  ray,  reflected  ray  and  transmitted  ray, 
respectively.  The  horizontal  components  of  B -field  in  Figure  (6a)  are  also  parallel  to  the 
boundary  and  therefore,  they  will  also  transmit  across  the  boundary  (consequence  of  Maxwell’s 
fourth  equation)  i.e., 


Bjcos(l)  —  Brcos(I)  =  Btcos(T ) 


(155) 


where  I  stands  for  angle  of  incidence,  which  is  same  as  angle  of  reflection.  T  stands  for  angle  of 
refraction.  Using  B  =  E/v  —  ( n/c)E ,  where  v  is  speed  of  light,  n  is  refractive  index  and  c  is  speed 
of  light  in  vacuum  in  Equation  (155)  results  in, 

—Ejcos(l) - -Ercos{I)  =  — Etcos(T )  (156) 

c  c  c 

where  n\  and  ni  represent  refractive  indices  of  incidence  medium  and  refraction  medium,  respec¬ 
tively. 

rtiEjCos(I)  — n\Ercos{I )  =  ri2Etcos(T )  (157) 

Equations  (154)  and  (157)  result  in, 


E,  n\cos(I)  —  n,2Cos[T) 
Ej  n\cos(I)  +  H2Cos(T) 


(158) 


Next  part  of  the  derivation  corresponds  to  Figure  (6b),  where  E-field  is  parallel  to  the  plane  of 
paper.  In  Figure  (6b),  B-field  is  perpendicular  to  the  plane  of  paper  and  parallel  to  the  boundary, 


(159) 


The  horizontal  component  of  E-field  is  also  parallel  to  the  boundary,  so  it  will  also  transmit  across 
the  boundary, 


Ejcos(l)  — Ercos(I )  =  Etcos(T ) 

(160) 

Using  the  relation,  B  =  ( n/c)E  in  eq.  (159), 

ri\  ri]  «2 

—Ej  -\ — -Er  =  —E, 
c  c  c 

(161) 

n\Ej  +  n\Er  =  n2Et 

(162) 

Equations  (160)  and  (162)  result  in, 

E,  n,2Cos(I) —  n\cos{T) 

Ej  n\cos{T)  +ri2Cos(I) 

(163) 

If  light  is  assumed  to  be  unpolarized,  the  reflectivity  or  reflectance  is  given  as, 


Substituting  Equations  (158),  (163)  into  Equation  (164)  and  using  normal  incidence, 
i.e.,  I  =  T  —  0,  /„  \  2 


R  = 


n2  —  n\ 
n2  +  rn 


(165) 
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Fresnel’s  equation  for  dielectric-conductor  interface  is  analogous  to  dielectric-dielectric  case.  772 
in  Eqs.  (158),  (163)  needs  to  be  replaced  with  complex  refractive  index  (772  —  ik2),  where  kn  is  an 
optical  parameter.  Fresnel’s  equation  for  dielectric-dielectric  interface  is  given  as  [37], 


Er  n\cos(I )  —  (772  —  ik2)cos[T ) 

Ej  n\cos(I)  +  (722  +  ik2)cos[T) 


(166) 


Er 

Ei 


(772  —  ik2)cos(I)  — n\cos(T ) 
77 1  co5(7’)  +  (772  —  ik2)cos(I ) 


(167) 


Based  on  Eqs.  (166),  (167),  the  s-polarized  and  p-polarized  reflectivities  are  given  as  [38], 


a2  +  b2  —2a  cos  I  +  cos2/ 

'  a2  +  b2  +  2a  cos  I  +  cos2 1 

a2  +  A>2  —  2a  sin/ tan/ +  sin2 /tan2/ 
s  a2  +  b2  +  2a  sin  /  tan  /  +  sin2 /  tan2 / 


a2  =  ^2  ^  (77^  —  k\  —  772  sin2  7)2  +  4t?^2  +  772  —  A:2  —  772  sin2 

b2  =  ^2  ^  (77^  —  k2  —  7? 2  sin2  7)2  +  4t7^2  —  772  -(-  A:2  +  77J  sin2/^ 

When  angle  of  angle  of  incidence  is  zero,  both  Rs  and  Rp  take  the  following  form, 

D  _  D  ^  +  (77l-  772)2 
^  P  k2  (77 !  772)2 

Furthermore,  with  the  assumption  of  sunlight  being  unpolarized,  reflectance  is  given  as, 

R  _  Rs  +  Rp  _  kl  +  (n\  772)2 

2  kj  +  (ni+n2)2 


(168) 

(169) 

(170) 

(171) 

(172) 

(173) 


Fresnel’s  Eq.  (173)  requires  values  of  refractive  indices  721,  n2  and  value  of  optical  parameter 
k2.  It  is  safe  to  assume  77 1  =  1,  which  corresponds  to  vacuum.  The  objective  is  to  find  refractive 
index  of  the  second  medium  n2  as  function  of  light  wavelength  or  frequency. 

A  case  of  an  electron  bound  to  its  nucleus  is  considered  and  it  is  trying  to  move  under  the  influence 
of  an  applied  electric  field,  E  [39].  An  electric  dipole  moment  is  created  because  of  electric  field’s 
attempt  to  separate  the  electron  from  its  nucleus.  The  system  can  be  described  as, 

mx  —  eE  —  kx  —  mjx  (174) 

where  m  is  mass  of  electron,  e  is  elementary  charge  and  k  is  spring  constant.  It  is  assumed  that 
electric  field  is  applied  in  x-direction,  ‘kx’  models  the  restoring  force  arising  from  attraction  of  nu¬ 
cleus,  mjx  is  sort  of  a  frictional  force  and  it  arises  from  collisions  of  electrons  with  other  particles. 

Resonance  frequency  is  given  as  (Dq  —  Eq.  (174)  can  be  written  as, 

mx  —  eE  —  m(OqX  —  myx  (175) 
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(176) 


x  +  yx  +  C0qX  —  — K 


m 


Cdo  =  0  corresponds  to  the  case  of  unbound  electrons  i.e.,  good  conductors  and  C0q  /  0  corresponds 
to  dielectrics,  y^  0  for  both  conductors  and  dielectrics.  The  electric  field  can  have  any  functional 
form,  and  it  will  be  assumed  of  the  form  E(t)  =  Ee2(0t.  Equation  (176)  has  the  standard  solution 
x(t)  =  xe2COt ,  where  the  amplitude  x  satisfies  the  relation, 


9  9  € 

-co x  +  jcoyx+conx  =  —  E 

m 


x 


'  / 

m 


cOq  —  co2  +  jcoy 


(177) 


(178) 


The  polarization  P  is  given  as, 


P  =  fVp  =  Ae.r  = 


—E 

m _ 

(On  —  co2  T  jcoy 


=  £oZ{co)E 


(179) 


£oZ(w)  - 


Nez 

m 


Oli 


co2  +  jcoy 


(180) 


where  N  is  the  total  number  of  elementary  dipoles  per  unit  volume,  Eq  is  permittivity  of  vacuum 
and  %(co)  is  electric  susceptibility  of  medium.  Now,  the  electric  flux  density  is  given  as, 


D  =  e(co)E  =  EqE  +  P 

e(co)E  =  e0(1  +  x(co))E 

e(co)  =  e0(l +*(«)) 


Nel 

m 


£{co)  =  £0H - Y  2  ,  • 

cOq-co-  +  jcoy 


(181) 

(182) 

(183) 

(184) 


where  e(co)  is  effective  permittivity.  Using  the  definition  of  ‘plasma  frequency’  as  co2  =  ^  in 
Equation  (184), 


,  ,  ,  eofflp 

e{co)  =  eQ  +  — 

cox  —  co2  T  jcoy 


(185) 


e(co)  will  have  a  real  part  and  an  imaginary  part,  so  let  e(co)  =  £'((o)  —  j£n(co).  e' (co)  is  char¬ 
acteristic  of  refractive  properties,  whereas  e!'(co )  is  characteristic  of  absorptive  properties.  On 
separating  Equation  (185)  into  real  and  imaginary  parts, 


E0COI(COI-C02) 

e((o)  =  £q  +  1 


(co2  -  co^  +  yco 2 


(186) 


e"(co)  = 


e  cot;  coy 


(co2  -  (Oq)  +  y2co2 


(187) 
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Around  the  resonance  frequency  top,  £r  shows  anomalous  behaviour  and  absorption  increases  sub¬ 
stantially.  A  small  detour  will  be  taken  here  to  have  a  look  at  the  case  of  dielectric  materials  as 
well.  Although  the  focus  of  this  research  is  conductors,  but  the  development  can  be  used  later  if 
dielectric  materials  are  ever  taken  up  for  research. 


In  actual  dielectric  materials,  there  exists  many  resonance  frequencies  corresponding  to  various 
vibrational  modes  and  polarization  mechanism  and  Equation  (185)  has  the  form. 


£(»>)  =  tO  +  3)E  fe?/('"’ei)) 


cof  -  to2  +  jay 

More  precise  application  of  quantum  mechanical  models  lead  to  essentially  the  same  equation 


(188) 


+  A0)g2/('»£o) 

(0Z  +  J(0Jji 


J>i  ^ 


(189) 


where  (Ojj  represents  transition  frequency  between  energy  levels,  Nj  and  Nj  correspond  to  popula¬ 
tion  of  lower  and  upper  energy  levels,  respectively,  fjj  are  called  oscillator  strengths.  On  using  the 


definition  of  refractive  index  n(co)  =  y  Equation  (188) 

2/  ,  ,  vn  Bi(of 

/72(m)  =  l+£- 


£ Of  -  ft)2  +  jCOJi 


(190) 


where  Bj  are  constants.  Usually,  the  above  equation  is  applied  to  regions  far  away  from  resonance 
so  that  Yi  =  0, 


7/  ,  ,  vn  Bj(0f 

R2  ft)  =  l+£  2  ' 

ft)“  -  ft)- 

(191) 

BjX2 

n  (“>=1+IA2_A2 

(192) 

This  is  known  as  Sellmeier  equation  and  is  valid  for  dielectric  materials.  The  values  of  the  wave¬ 
lengths  corresponding  to  resonance  condition  and  the  constants  Bj  are  obtained  experimentally. 


For  conductors,  (Oq  =  0,  and  thus  Equation  (185)  can  be  written  as, 


e(co) 

£o 


=  1  + 


ft); 


j(o{y+j(o) 


(193) 


The  above  model  is  referred  to  as  Drude  Model.  A  more  detailed  modelling  using  interband  terms 
gives  rise  to  Lorentz-Drude  model,  according  to  which  the  complex  dielectric  function  is  given  as 

[41], 


£,-{(0)  =  =  £rf)(C0)  +  £rh\<o) 


£o 


(194) 


where  £^\(0)  corresponds  to  free-electron/Drude  model  and  £,h>  { ft))  corresponds  to  interband 
part  or  Lorentz  model  i.e., 


£rf\co)  =  1  - 


Q.2 

P 


co(co-iTo) 


(195) 
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(196) 


>(b) 


(®)  =  £ 


£  K 


fj% 

co2)  +  icoFj 


where  (Op  is  plasma  frequency,  1  is  the  number  of  oscillators  with  frequency  C0j,  strength  fj,  and 
lifetime  1  /Fj,  and  Flp  —  \/JoCQp  is  plasma  frequency  associated  with  intraband  transitions  with 
oscillator  strength  /o  and  damping  constant  Tq. 


A  modified  version  of  classical  Lorentz-Drude  model  is  proposed  by  Brendel  and  Bormann  [40] 
according  to  which  the  complex  dielectric  function  is  given  as. 


£r((0)  =  1 


Qt 


co(tO-rTo) 


/ 


+  £Xj((o) 

7=1 


(197) 


Xj(0)) 


1 

\/2nOj 


M _ dx 

co2)  +  icoFj 


(198) 


where  (Op  =  ^  is  plasma  frequency,  with  N,  e,  £o,  m  being  number  of  conduction  electrons  per 

unit  volume,  electronic  charge,  permittivity  of  vacuum  and  electronic  mass,  respectively.  1  is  the 
number  of  oscillators  with  frequency  C0j,  strength  fj,  and  lifetime  1  /Fj,  and  Flp  =  \ff)(Op  is  plasma 
frequency  associated  with  intraband  transitions  with  oscillator  strength  /o  and  damping  constant 
To-  <Jj  is  standard  deviation  associated  with  the  Gaussian  modelling  of  the  dielectric  function  and 
its  value  is  experimentally  obtained,  co  is  angular  frequency  of  light. 


Rakic  [41]  has  given  experimental  tables  for  Brendel-Bormann  model  parameters  /o,  fj.  To,  Fj, 
C0j,  <Jj  for  some  commonly  used  materials  including  the  material  Aluminum,  which  is  used  in  this 
research  work. 


Analytical  solution  for  XjiP)  exists  and  is  given  as, 


Xj  = 


ifjtoj 

I'JlcijOj 


U 


1/2, 1/2,- 


aj  ~  °fj 

V20; 


2i 


+  u 


2-i 


(199) 


Here,  U  is  known  as  Kummer  function  of  second  kind  and  17(1/2, 1/2,  z2)  —  \fjtc:r  erj'c(z),  where 
erfc(x)  =  ff  e~rdt  is  the  complementary  error  function.  The  parameter  aj  is  given  by  its  real 
and  imaginary  parts, 


where, 


cij  =  a'j  +  ia" 

4  =  ^((1+(r,»2)1/2+1)1/" 
4'  =  ^((1  +  (rJ/“)2)1/2-1)1/2 


(200) 

(201) 

(202) 
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Now,  the  complex  refractive  index  is  given  as, 


h  —  n  —  jk 


(203) 


where  n  is  the  usual  real  refractive  index,  k  is  extinction  coefficient.  Now,  the  square  of  amplitude 
of  this  complex  refractive  index  ( n )2  —  (n2  —  k2)  —  2  jkn  is  equal  to  the  complex  dielectric  function. 

Thus,  if  complex  dielectric  function  £,-(&))  —  en  —  j£r2,  then  from  £r(co)  —  (n)2, 

£n  =n2  —  k2  (204) 

£n  —  2kn  (205) 


Solving  Equations  (204)  and  (205), 

n  -  ^  \/£n  +  \fen  +£n 

k=^  \J~£ri  +  \/en  +  e£ 


(206) 

(207) 


3.9  Capacitance  Modeling 

In  the  previous  section,  currents  for  determining  voltage,  are  modeled.  In  order  to  compute  object 
charge  (to  be  subsequently  used  for  computing  Lorentz  force)  from  voltage,  Equation  (6)  requires 
body  capacitance,  whose  modeling  depends  upon  the  orbit  regime.  Primary  reason  for 
difference  in  capacitance  modeling  arises  from  difference  in  Debye  length  relative  to  body 
dimension.  For  an  object  with  dimension  ranging  from  tens  of  centimeters  to  few  meters,  Debye 
length  is  much  larger  than  body  dimension  in  GEO  region,  whereas  Debye  length  is  much  smaller 
than  body  dimension  in  LEO  region.  As  a  consequence  to  this,  certain  first  order 
approximations  can  be  made  in  the  capacitance  formulation,  as  will  be  seen  shortly. 


3.10  Spherical  Capacitance  Model 

For  two  concentric  conducting  spheres  [28]  having  radii  a  and  b  respectively, 

Q 


AV  — 


b  1 

^)dR 


4/rt'o  Ja  R2 


AV  = 


Q  f  1  1 

47T£o  \a  b 


Thus,  the  capacitance  C  is, 


C  = 


4^£0 
T  P 

a 


ri- 1 

\a  b 


(208) 

(209) 

(210) 


where  £o  is  vacuum  permittivity.  At  GEO  altitude,  a  =  R  and  b  =  L,  where  R  is  radius  of  the 
assumed  spherical  conducting  debris  and  L  is  the  Debye  length. 


C  = 


AtceqR 


(211) 
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With  j-  «  1  for  GEO,  Equation  (211)  can  be  approximated  as, 


C  =  4  7l£0Rz  (7  +  7 

,  K  L 


(212) 


Chunshi  [42]  makes  further  assumption  by  neglecting  I  compared  to  i,  and  Equation  (212)  becomes, 


R  =  4  tieqR 

Using  a  =  R  and  b  =  (R+L)  for  LEO,  Equation  (21 1)  results  in, 

c_  47T£q R{R  +  L) 


(213) 


(214) 


With  j-  »  1  for  LEO,  Equation  (214)  can  be  approximated  as, 


C  = 


4  keqR2 


(215) 


The  Debye  length  L  can  be  obtained  from  the  following  equation, 

1  y(nJq2j\ 

L2  L\£okTj) 


(216) 


where  n,  q,  T  and  k  denote  number  density  of  plasma  particles,  charge  of  particle,  characteristic 
temperature  of  particle  distribution  and  Boltzmann  constant,  respectively,  for  particle  type  j. 


3.11  Analytic  Expressions  for  Orbital  Perturbation  Due  to  Lorentz  Force 

Numerical  simulation  is  time  expensive  and  does  not  explicitly  highlight  the  factors  effecting  the 
perturbations  in  various  orbital  elements.  This  section  is  dedicated  to  developing  analytic  expres¬ 
sions  for  perturbations  in  orbital  elements  over  one  orbital  period  time  using  Gauss’s  variational 
equations.  The  expressions  developed  apply  to  a  Keplerian  orbit  of  small  eccentricity  perturbed  by 
Lorentz  force. 


The  classical  Gauss’s  variational  equations  representing  time  rate  of  change  of  orbital  elements, 
are  given  as  [14], 


da 

dt 

de 

dt 


2 

n\J  1  —  e2 

Vl  —  e2 
na 


Se  sin  /  +  T  ( 1  +  e  cos  /) 
S  sin  /  +  T  (cos  /  +  cos  E) 


di  r  cos  (co  +  f) 

—  = -  J  W 

dt  na2V  1  —  e2 

dD.  r  sin  (co  +  f)  w 

dt  na2V  1  —  e2s\ni 


(217a) 

(217b) 

(217c) 

(217d) 


35 

Approved  for  public  release;  distribution  is  unlimited. 


—  Scos/  +  t(i  +  —  ^  sin/  —  cosz^ 


(217e) 


where  a,  e,  i,  fi.  ft),  f  represent  the  Keplerian  orbital  parameters  of  semi-major  axis,  eccentricity, 
inclination,  right  ascension  of  ascending  node  (RAAN),  argument  of  perigee  and  true  anomaly, 
respectively,  n,  r,  p,  E  represent  mean  motion,  radial  distance,  semilatus  rectum  and  eccentric 
anomaly,  respectively.  S  and  W  are  the  components  of  Lorentz  acceleration  in  radial  direction  and 
perpendicular  to  orbital  plane  direction,  respectively.  T  represents  component  of  Lorentz  acceler¬ 
ation  in  a  direction  that  completes  the  right  handed  orthogonal  system.  Peng  et  al.  [14]  calls  S, 
T,  W  frame  as  local-vertical-local-horizontal  (LVLH)  frame.  Obtaining  expressions  for  S,  T,  W  in 
terms  of  orbital  parameters  a,  e,  i,  Q,  ft),  f,  as  given  in  Peng  et  al.  [14],  is  the  key  to  integrating 
Equations  (217a)-(217e). 


According  to  classical  Electrodynamics,  Lorentz  acceleration  of  an  object  is  given  as, 


/  q 

cil=  — 
V  m 


—  {(be  xr)  xB  +  rxB 


(218) 


where  q  is  charge  of  the  object,  m  is  mass  of  the  object,  (be  is  the  angular  velocity  vector  of  Earth,  r 
is  the  geocentric  position  vector  of  the  object,  B  is  Earth’s  magnetic  field  vector  and  r  is  the  inertial 
velocity  vector  of  the  object. 

Use  of  the  relation  B  —  ^  [3  ( N  ■  f )  r  —  N]  and  Earth-centered  inertial  (ECI)  vectors  a)e  —  [0  0  (Oe]T , 
r—[xy  z]T,  N  =  [Nx  Ny  1VZ] r,  r—[xy  z]T ,  r=[f  y-  f\T  in  Eq.  (218),  where  Bq  is  magnetic  dipole 
moment  of  Earth,  r  is  the  geocentric  distance  of  the  object,  r  represents  unit  vector  in  the  radial 
direction  and  N  is  the  unit  vector  in  the  direction  of  dipole,  results  in  ECI  Lorentz  components 

[14], 


/  (j  \  /g  \  / 

dLx=  (  —  )  (  -^  )  [3x(zy -yz-xz(De)\Nx{t)  +  [3yzy - (3y2 - r2)z-3xyz(De]Ny(t)  +  [{3z2 - r2)y 
\  m  /  V  r  /  \ 


3 yzz  -  COex(3z2  -  r2)]Nz  )  (219) 


ftLv =  (  —  )  (  ~ r  )  (  [{3x2 - r2)z-3xzx-3xyz(oe]Nx{t)  +  [3y{xz- zx-yz(oe)\Nv{t)  +  [3xzz-  {3z2 
\m /  \  '  /  \ 


r2)x-(Oey{3z2-r2)]Nz\  (220) 
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Peng  et  al.  [14]  uses  two  coordinate  transformations:  the  first  transformation  is  performed  from 
ECI  coordinate  system  to  an  intermediate  frame  known  as  instantaneous  nodal  inertial  coordinate 
frame,  and  the  second  transformation  is  performed  from  this  intermediate  frame  to  local-vertical- 
local  horizontal  (LVLH)  frame.  The  instantaneous  nodal  inertial  coordinate  system  ( OX'Y'Z ')  is 
defined  as  follows:  X'  lies  along  instantaneous  intersection  of  equatorial  and  orbital  planes,  Z'  is 
perpendicular  to  equatorial  plane  and  Y'  completes  the  right-handed  system.  OX'Y'Z'  is  obtained 
by  rotating  the  ECI  frame  by  Q  about  Z-axis.  One  can  transform  OX'Y'Z'  frame  to  LVLH  frame 
by  performing  following  rotations  in  the  order:  angle  i  about  X'  axis,  angle  0  =  (ft)  +  /)  about 
newly  formed  Z'new  axis. 

Following  the  two  transformations,  Lorentz  acceleration  components  aix ,  aiy,  aiz  in  ECI  frame 
can  be  expressed  in  terms  of  Lorentz  acceleration  components  S,  T,  W  in  LVLH  frame  as  [14], 
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a  j  v 


where  0  =  (ft)  +  /)  is  the  argument  of  latitude. 

The  dipole  direction  components  Nx,  Ny,  Nz  are  given  as, 


~n; 
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(223) 
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where  8m  represents  angle  between  dipole  north  pole  and  geographic  north  pole.  9m  is  sometimes 
referred  to  as  coelevation  of  the  dipole.  am  is  given  as, 

ccm  —  (Xqq  +  (Det  +  (j),n  —  £2  (224) 

where  OCco  represents  right  ascension  of  Greenwich  at  reference  time  (time  of  start  of  orbit  propa¬ 
gation  in  simulation)  and  <j>m  represents  the  east  longitude  of  dipole.  Bq,  dm  and  <j)m  can  be  obtained 
using  following  relations  from  Wertz  [16], 

fio  =  4\/(*?)2  +  («!)2  +  ('>!)2 

8°l 

cos  G,„  =  , 

\!  {gy- + {g\)2 + (h\y~ 

,  h\ 

tan  <j)m  =  -[ 

Si 

where  g®,  gj,  h\  represent  IGRF  coefficients  and  Re  is  mean  radius  of  Earth 

Equations  (222)  -  (227)  can  now  be  used  in  Equation  (217),  which  can  then  be  integrated  to  obtain 
one-period  orbital  perturbations.  Peng  et  al.  [14]  derived  expressions  for  one-period  orbital  per¬ 
turbations,  but  they  assumed  that  the  inertial  direction  of  Earth’s  magnetic  dipole,  Nx(t)  and  Ny(t), 
remain  constant  during  one-orbital  period.  This  assumption  is  sufficiently  accurate  for  a  low  Earth 
orbit  whose  period  is  much  small  compared  to  Earth’s  rotation  period.  The  formulas  derived  in 
this  report  for  one-orbital  period  perturbations  work  well  for  all  altitudes.  In  this  research,  fol¬ 
lowing  approximations  are  made  while  deriving  the  perturbation  expressions:  (1)  eccentricity  is 
small  (near- zero),  which  allows  for  approximation  of  time  as  t  ~  „  (2)  charge-to-mass  ratio  of  the 
body  is  constant  (3)  the  factor  (1  +  ecos /)-1,  which  appears  in  the  derivation  is  approximated  by 
second-order  Taylor  series  (4)  the  orbital  elements  [a,  e,  i,  ft),  £2]  are  assumed  to  remain  constant 
over  the  integration  period. 


(225) 

(226) 

(227) 


In  obtaining  the  one-orbital  period  perturbations,  true  anomaly  is  considered  as  the  free  variable, 


A<=f 


Jo  V  dt  df 


(228) 


where  x  =  [a,  e,  i,  ft),  £2]. 


Three  classes  of  perturbation  formulas  are  developed, 

1.  General  formulas:  applicable  for  all  orbit  regions  (except  the  precise  GEO  altitude  with 
c oe  —  n,  because  of  the  factor  (^  —  1)  in  denominator). 

2.  LEO  forumlas:  applicable  only  for  LEO  orbits,  and  are  derived  from  general  formulas  using 
the  approximations  (Oe  «n  and  e2  «  1. 

3.  GEO  formulas:  applicable  for  GEO  orbit  and  orbits  around  GEO  altitude,  and  are  derived 
from  general  formulas  taking  the  limit  %  — >  1  ,  where  %  — 


38 

Approved  for  public  release;  distribution  is  unlimited. 


Following  relations  are  useful  in  the  derivation, 


p  —  a(l  —  e2) 

(229a) 

"  =  V? 

(229b) 

e  A  cos  f 

cos  F  — 

lAecos  / 
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r  = - 

1  Aecos  / 

(229d) 

dt  1  /  r\2  1 

df  n\a)  yj 1-e 2 

(229e) 

3.12  General  and  LEO  Perturbation  Expressions  for  One  Orbital  Period 

The  expressions  for  general  perturbation  formulas  and  LEO  perturbation  formulas  have  identical 
form,  but  different  coefficients  need  to  be  inserted.  The  coefficients  for  general  perturbation  for¬ 
mulas  can  be  obtained  from  Appendix  A,  and  the  coefficients  for  LEO  perturbation  formulas  can 
be  obtained  from  Appendix  B. 
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(230) 


In  comparison  with  the  formulation  of  Peng  et  al  [14]  that  can  be  found  Appendix  C,  the  non¬ 
constant  dipole  leads  to  further  terms,  in  the  power  of  the  eccentricity.  The  simple  sini  dependence 
is  extended  with  the  varying  dipole  direction. 
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Comparable  as  in  the  case  of  the  semi-major  axis,  relaxation  of  the  eccentricity  and  dipole  stability 
leads  to  higher  order  terms  in  the  eccentricity  and  the  inclination  dependence  is  contrasted  with  the 
dipole  orientation  terms. 
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2  k  ox, 
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In  comparison  to  Peng,  higher  order  terms  in  the  eccentricity  occur.  The  J  terms  are  not  referring 
to  gravitational  constants  but  can  be  found  in  the  appendix  as  well. 
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Again,  higher  order  terms  do  occur.  The  term  labelled  N  are  placeholders  and  can  be  found  in  the 
appendix.  They  are  unconnected  to  the  dipole  direction. 
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3.13  GEO  Perturbation  Expressions  for  One  Orbital  Period 


In  the  following,  the  orbital  averaged  expressions  of  the  Lorentz  force  induced  perturbations  of 
each  of  the  orbital  elements  in  the  geostationary  region  are  listed. 
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Compared  to  Peng  et  al  [14],  as  in  the  previous  case,  higher  order  terms  are  picked  up  in  eccentric¬ 
ity.  The  dipole  dependence  of  the  Lorentz  force  in  combination  with  the  time  dependence  of  the 
dipole  leads  to  a  coupled  term  of  the  sini  and  the  pole  change. 
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In  the  eccentricity,  the  overall  structure  of  the  terms  is  comparable  to  Peng,  however,  the  time 
varying  dipole  terms  occur. 
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In  general  the  terms  for  the  geosynchronous  region  are  simpler  and  easier  to  control  compared 
to  the  low  earth  orbital  region  or  the  more  general  expressions.  This  is  due  to  the  fact  that  the 
approximation  limit  of  —  could  be  utilized. 


4.0  RESULTS  AND  DISCUSSION 


Several  orbits  corresponding  to  LEO  and  GEO  regions  are  simulated.  For  all  the  simulation  cases, 
an  aluminum  sphere  of  1  nr  surface  area  is  considered.  Simulation  time  is  four  days.  Epoch 
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of  initialization  is  12:00:00  UTC,  March  21,  2000.  Earth’s  gravity,  Sun  and  Moon  gravitational 
perturbation,  direct  solar  radiation  pressure  (SRP),  atmospheric  drag  (only  for  LEO)  and  Lorentz 
force  perturbation  are  included  in  the  modeling.  Gravitational  potential  has  been  modeled  up  to 
degree  and  order  12.  Atmospheric  drag  coefficient  is  taken  as  two  for  spherical  body.  SRP  diffuse 
reflection  coefficient  is  taken  as  .035.  Dipole  model  has  been  chosen  for  Earth’s  magnetic  field. 
Earth’s  shadow  is  modeled  as  a  cylinder.  For  LEO  simulations,  plasma  electrons  are  assumed  to 
be  collected  over  entire  surface  area,  whereas  photons  and  plasma  ions  are  assumed  to  be  collected 
over  frontal  surface  area.  For  GEO  simulations,  plasma  electrons  and  ions  are  assumed  to  be  col¬ 
lected  over  entire  surface  area,  whereas  photons  are  assumed  to  be  collected  over  frontal  surface 
area.  This  difference  between  GEO  and  LEO,  as  stated  earlier,  is  because  LEO  ions  are  at  rest 

relative  to  body  and  hence  ram  into  frontal  surface  area.  Two  area-to-mass  ratio  objects  are  inves- 

2 

tigated,  one  corresponds  to  general  satellites  with  low  area-to-mass  ratio  (LAMR)  of  .02  and 

2 

the  other  corresponds  to  high  area-to-mass  ratio  (HAMR)  of  23.6  'j^. 

4.1  Low  Earth  Orbit  Simulation 

Three  different  groups  of  LEO  are  simulated.  The  first  LEO  group  is  characterized  by  low  initial 
inclination;  the  second  LEO  group  is  characterized  by  high  initial  inclination;  and  the  third  LEO 
group  comprises  of  sun-synchronous  orbits.  Table  3  describes  the  LEO  cases  that  are  simulated. 

Table  3.  The  LEO  Cases 


a  =  6778  Km.  ftj  =  30°.  fl  =  60°,  v  =  100°.  e  =  .001 


LEO  group  1  (i  =  .  1°)  LEO  group  2  (i  =  89°)  LEO  group  3  (i  =  97.0292°) 


AMR  =  .02 

AMR  =  .02 

AMR  =  .02 

AMR  =  23.6 

AMR  =  23.6 

AMR  =  23.6 

Two  plasma  conditions  are  simulated,  one  resulting  in  low  body  charge  and  the  other  leading  to 
high  body  charge.  For  low-charge  plasma  condition,  International  Reference  Ionosphere  (IRI)  is 
used  to  obtain  electron  temperature,  ion  temperature,  electron  density,  percentage  ion  composition 
for  a  given  date,  latitude,  longitude  and  height  profile.  Composition  of  0+,  H+,  He+,  Ot,  NO+ 
and  N+  ions  can  be  obtained.  High-charge  plasma  condition  corresponds  to  body  entering  auroral 
ovals  in  the  north  and  south  geomagnetic  hemispheres.  Location  of  auroral  oval  is  determined 
using  Holzworth’s  interpolation  curve  [43], 

0  =  A]  +  A2  cos  (0  +A3)  +A4COS  (2<j>  +2A5)  +A6cos(30  +3A7)  (239) 

where  0  is  corrected  geomagnetic  co-latitude  and  (j)  represents  magnetic  local  time  in  radians.  A 
table  listing  values  of  coefficients  A\  through  A7  for  poleward  edge  and  equatorward  edge  of  the 
auroral  oval  is  given  by  Holzworth  [43]  at  various  geomagnetic  conditions  discriminated  by  the 
parameter  Q  (IGY  quarter-hourly  range  index  of  magnetic  activity)  values.  For  the  simulations  in 
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this  report,  Q  is  two  (relatively  quiet  geomagnetic  activity).  Corrected  geomagnetic  co-latitude  and 
magnetic  local  time  are  obtained  using  a  software  package  developed  by  Shepherd  [44].  The  high- 
charge  plasma  densities  and  temperatures  inside  auroral  ovals  are  simulated  as  Spenvis  Fontheim 
with  ECSS  cold  background  environment  [45],  which  is  defined  in  Table  4. 


Table  4.  Plasma  Parameters  Inside  Auroral  Oval 


Particles 

Number  Density  {cm  3) 

Energy  (eV) 

0+  ions 

125 

.2 

electrons  1 

125 

.2 

electrons2 

1.482 

12940 

Figure  7  shows  the  voltage  and  current  plots  of  LEO  LAMR  bodies  Under  Low-charge  plasma 
condition.  Figure  7a  plots  the  voltage  comparison  of  group  1,  group  2,  group  3  LEO  objects.  The 
body  voltage  plots  have  a  short-period  variation  of  one  orbital  period  and  a  long-period  variation  of 
one  day.  These  periodic  variations  in  voltage  result  from  periodicity  in  altitude  or  the  LEO  plasma 
data  dependence  on  altitude;  the  altitude  variations  arise  mainly  because  of  the  orbit  eccentricity 
value  and  its  variation.  Figure  7b  shows  the  logarithm  of  absolute  value  of  current  for  0.1°  inclined 
LEO.  Plasma  electron  current  is  the  dominating  current,  and  hence  the  equilibrium  voltages  are 
negative.  In  the  photoelectric  current  plot  of  Figure  7b,  when  an  object  is  in  Earth’s  shadow, 
the  photoelectric  current  is  zero  (or  its  logarithm  value  is  not  defined),  whereas  it  is  a  constant 
negative  value  when  the  object  is  in  sunlight.  Since  photoelectric  current  is  several  orders  of 
magnitude  smaller  than  other  currents,  the  voltage  changes  resulting  from  sunlight-to-shadow  and 
shadow-to- sunlight  transitions  are  not  significant.  Figure  7c  shows  the  difference  between 
currents  of  group  2  and  group  1  (for  truncated  duration:  orbit  30  -  orbit  32),  Figure  7d  shows 
difference  between  currents  of  group  3  and  group  1.  As  can  be  seen,  group  2  and  group  3 
currents  are  of  similar  order  as  that  of  currents  of  group  1.  Difference  of  photoelectric  currents  is 
deliberately  not  included  in  Figure  7c  and  Figure  7d,  as  different-group  orbits  enter  or  leave 
Earth’s  shadow  region  at  different  times;  the  values  of  photoelectric  currents  for  group  2  or  group 
3,  are,  however,  same  as  that  of  group  1. 
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Figure  7.  Charging  Characteristic  of  LEO  Under  Low-Charge  Plasma  Conditions  (LAMR) 


Figure  8  shows  the  voltage  for  LEO  objects  Under  High-charge  plasma  conditions.  As  compared 
to  low-charge  plasma  condition  depicted  in  Figure  7,  where  the  maximum  negative  voltage  is 
less  than  unity,  here  the  voltage  can  reach  thousands  of  negative  volts.  The  body  voltages  inside 
auroral  oval  reach  approximately  -2.6K  Volts  under  sunlit  condition  and  reach  approximately 
-60K  Volts  (orbit  number  51  and  53  in  Figure  8a,  e.g.),  when  in  Earth’s  shadow.  When  outside 
the  auroral  oval,  body  voltages  are  small  and  of  comparable  order  as  that  of  low-charge  plasma 
conditions. 

Figure  9  shows  the  common  logarithm  of  magnitudes  of  perturbing  accelerations  for  FEO  orbits. 
The  most  dominant  and  the  least  dominant  perturbing  accelerations  correspond  to  higher  harmon¬ 
ics  of  Earth’s  gravitational  force  and  Forentz  force,  respectively.  The  jumps  in  Forentz  force  in 
higher  inclined  orbits  are  due  to  B  nearly  aligning  up  with  ~vrei. 
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(b)  Body  Voltage  (Group  3) 


Figure  8.  Voltage  Characteristic  of  LEO  Objects  (LAMR,  High-Charge  Plasma 

Condition) 


In  order  to  investigate  the  importance  of  Lorentz  force,  difference  between  simulation  cases 
including  Lorentz  force  and  simulation  cases  excluding  Lorentz  force  for  low-charge  plasma 
conditions  are  displayed  in  Figures  10,  11,  12.  Figure  10  shows  the  in-track  perturbations.  For 
the  in-track  perturbation,  the  dynamics  is  clearly  secular  for  all  the  three  groups.  In  addition  to 
the  long-term  secular  variation,  there  is  also  a  short-period  variation,  most  clearly  visible  in 
group  1  in-track  plots;  the  periodic  variation  corresponds  to  the  body  voltage  periodic  variation 
as  seen  in  Figure  7a.  The  positive  values  of  in-track  perturbations  also  suggest  that  the  position 
in  Lorentz  force  included  orbit  is  always  ahead  of  the  position  in  Lorentz  force  excluded  orbit. 
The  in-track  perturbations  are  of  the  order  of  10-3  m  for  all  the  three  groups. 

Figure  1 1  shows  the  cross-track  perturbations.  For  cross-track  perturbations,  the  mean  is  approxi¬ 
mately  zero,  which  implies  that,  on  an  average,  there  is  no  change  in  orbital  plane  due  to  Lorentz 
force.  For  all  the  three  groups,  there  is  a  short-period  variation  with  time-scale  of  one  orbital  pe¬ 
riod  and  increasing  amplitude.  A  factor  contributing  to  the  increasing  short-period  amplitude  of 
cross-track  perturbations  is  the  increasing  inclination  of  Lorentz  included  orbit  relative  to  Lorentz 
excluded  orbit  over  time.  For  higher  inclined  group  2  or  group  3  orbits,  there  is  also  a  long-period 
variation  with  time-scale  of  one  day.  Cross-track  perturbations  are  also  influenced  by  orbit  incli¬ 
nation,  as  the  perturbations  are  an  order  higher  for  higher  inclined  group  2  or  group  3  orbits  as 
compared  to  lowly  inclined  group  1  orbit,  whose  cross-track  perturbations  are  of  the  order  of  10  4 
to  10-3  m. 


Figure  12  shows  the  radial  perturbations.  For  radial  perturbations  of  lowly  inclined  group  1  orbit, 
the  mean  is  approximately  a  negative  constant  value.  Radial  perturbations  of  group  1  orbit  exhibit 
a  short-period  variation  of  one  orbital  time  period  with  an  increasing  amplitude.  Possible  explana¬ 
tions  for  the  increasing  amplitude  include  the  increasing  amplitude  of  variation  in  the  semi-major 
axis  difference  between  Lorentz  included  and  Lorentz  excluded  orbits.  Higher  inclined  group  2  and 
group  3  orbits  exhibit,  on  an  average,  a  decreasing  secular  behavior.  A  factor  contributing  to  this 


46 

Approved  for  public  release;  distribution  is  unlimited. 


secular  behavior  in  radial  perturbations  is  the  decreasing  semi-major  axis  of  Lorentz  force  included 
orbit  relative  to  semi-major  axis  of  Lorentz  force  excluded  orbit  over  time.  A  short  one-orbital  pe¬ 
riod  variation  and  a  long  one-day  period  variation  are  also  exhibited  by  the  radial  perturbations  of 
group  2  and  group  3  orbits.  The  radial  perturbations  are  of  the  order  of  1 0  5  m  for  group  2  and  d3 
and  of  10  4  m  for  group  1. 


0  10  20  30  40  50  60 

Number  of  Orbits 


0  0.5  1  1.5  2  2.5  3  3.5  4 

Number  of  Days 

(a)  LEO  Group  1 


Number  of  Days 

(c)  LEO  Group  3 


0 


-14 - ' - ' - ' - 1 - 1 - »- 

0  10  20  30  40  50  60 

Number  of  Orbits 


0  0.5  1  1.5  2  2.5  3  3.5  4 

Number  of  Days 

(b)  LEO  Group  2 


- Sun  3rd  Body  Perturbation  LAMR 

- Moon  3rd  Body  Perturbation  LAMR 

Solar  Radiation  Pressure  LAMR 

- Lorentz  Force  LAMR 

Higher  Harmonics  of  Earth  Gravity  LAMR 
— Atmospheric  Drag  LAMR 


Figure  9.  Common  Logarithm  of  Perturbing  Acceleration  Magnitudes  for  LEO 
Orbits  Under  Low-Charge  Plasma  Conditions 
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Figure  10.  In-track  Perturbations  for  Lorentz  Force  Perturbed  Low  Earth  Orbit  Relative 
to  Lorentz  Force  Excluded  Orbit  Under  Low-Charge  Plasma  Conditions  (LAMR  Object) 
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Figure  11.  Cross-track  Perturbations  for  Lorentz  Force  Perturbed  Low  Earth  Orbit 
Relative  to  Lorentz  Force  Excluded  Orbit  Under  Low- Charge  Plasma  Conditions 

(LAMR  Object) 
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Figure  12.  Radial  Perturbations  for  Lorentz  Force  Perturbed  Low  Earth  Orbit  Relative 
to  Lorentz  Force  Excluded  Orbit  Under  Low-Charge  Plasma  Conditions  (LAMR 

Object) 

Figure  13  shows  in-track,  cross-track  and  radial  perturbations  between  Lorentz-included  and 
Lorentz-excluded  simulation  cases  for  LEO  LAMR  objects  Under  High-charge  plasma 
condition.  High-charge  plasma  conditions  are  created  through  auroral  passages,  which  do 
not  apply  to  lowly  in-clined  group  1  orbit.  Compared  to  low-charge  plasma  conditions  given 
by  Figures  10,  11,  12,  the  perturbations  in  Figure  13  are  several  order  magnitudes  larger  due 
to  high  negative  voltages  inside  auroral  oval  as  seen  in  Figure  8.  The  in-track  and  cross-track 
perturbations  are  of  the  order  of  two  to  3.5  meters  for  group  2.  Group  3  experiences  a  large 
secular  trend  in  the  in-track  direction  of  up  to  35  meters  and  a  periodic  perturbation  of  up  to 
5  meters  in  the  cross  track  direction.  The  radial  perturbations  are  of  the  order  of  0.25  m 
for  group  2  and  of  four  meters  for  group  3.  The  secu-lar/periodic  trend  is  similar  to  that  of 
low-charge  plasma  condition.  The  steep  change  in  group  2  in-track  and  radial  perturbations 
around  orbit  51  is  due  to  large  negative  body  voltage  as  observedin  Figure  8a,  same  holds  for  the 
step  function  like  increase  for  group  3  (compare  with  Figure  8b). 
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Figure  13.  In-Track  Perturbations,  Cross-Track  perturbations,  Radial  Perturbations  for 
Lorentz  Force  Perturbed  Low  Earth  Orbit  Relative  to  Lorentz  Force  Excluded  Orbit  Under 

High-Charge  Plasma  Conditions  (LAMR  object) 
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HAMR  LEO  objects  (initial  a  =  6778  Km)  crash  into  Earth  after  about  five  orbital  time-periods 
because  of  high  atmospheric  drag  forces.  Tables  5,  6  show  a  single  orbit  in-track,  cross-track 
and  radial  perturbations  for  HAMR  objects  in  LEO  Under  Low-charge  and  high-charge  plasma 
conditions  respectively. 


Table  5.  Perturbations  of  HAMR  Object  in  LEO  Under  Low-Charge  Plasma  Conditions 

After  One  Orbital  Time-Period 


Orbit  Group 

In-track  Perturbations  (m) 

Cross-track 

(m) 

Perturbations 

Radial  Perturbations  (m) 

Group  1 

1.39  x  KU1 

i  x  i<r2 

-7  x  1(T3 

Group  2 

2.7  x  ltr2 

—3.9  x  1(U2 

—4  x  ltr3 

Group  3 

6  x  ltr3 

—4  x  ltr2 

-3  x  ltr3 

Table  6.  Perturbations  of  HAMR  Object  in  LEO  Under  High-Charge  Plasma  Conditions 

After  One  Orbital  Time-Period 

Orbit  Group 

In-track  Perturbations  (m) 

Cross-track 

(m) 

Perturbations 

Radial  Perturbations  (m) 

Group  2 

5.34 

-1.91  x  101 

_4.49  x  1(T2 

Group  3 

5.85  x  KT1 

-1.63  x  101 

-3.04  x  10~2 

4.2  Geosynchronous  Earth  Orbit  Simulation 

Three  different  groups  of  geosynchronous  orbits  are  considered;  the  first  GEO  group  is  character¬ 
ized  by  low  inclination  (.1°)  and  small  eccentricity  (.001)  as  initial  parameters;  the  second  GEO 
group  is  characterized  by  inclination  of  15°  and  eccentricity  of  .015;  and  the  third  GEO  group  is 
characterized  by  inclination  of  40°  and  eccentricity  of  .015.  The  initial  argument  of  perigee  and 
right  ascension  of  ascending  node  values  are  kept  same  as  that  of  LEO  simulations.  The  initial  true 
anomaly  values  are  taken  such  that  the  initial  longitude  value  is  75°  E  (GEO  libration  point).  The 
GEO  cases  are  given  in  Table  7. 
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Table  7.  The  GEO  Cases 


a  =  42164  Km.  (0  =  30°.  G  =  60° 

GEO  group  1 

GEO  group  2 

GEO  group  3 

(i  =  .1°,  e  =  .001,  v  =  344.15°) 

(i=  15°,  e  =  . 015,  v  =  344.6) 

(i  =  40°.  e  =  .015,  v  =  349°) 

AMR  =  .02 

AMR  =  .02 

AMR  =  .02 

AMR  =  23.6 

AMR  =  23.6 

AMR  =  23.6 

Similar  to  LEO,  two  different  plasma  conditions  are  simulated,  one  corresponding  to  low-charge 
conditions  and  the  other  corresponding  to  high-charge  conditions.  Denton  [46]  provides  GEO 
plasma  parameters  for  three  classes  of  particles:  hot  electrons  (30  eV  to  45  KeV),  hot  protons  (100 
eV  to  45  KeV)  and  low  protons  (1  eV  to  100  eV).  Data  corresponding  to  one-solar-cycle  (1990  - 
2001)  averaged  low  proton  density,  hot  proton  density,  hot  electron  density,  hot  proton  perpendic¬ 
ular  and  parallel  temperatures,  hot  electron  perpendicular  and  parallel  temperatures  (Maxwellian 
temperatures  for  distributions  in  the  perpendicular  and  parallel  directions  of  magnetic  field,  respec¬ 
tively)  can  be  obtained  from  Denton  [46] .  The  net  temperature  for  the  overall  distribution  is  then 
taken  as 

T  —  +  ^-7j|  (240) 

Temperature  distribution  for  low  protons  is  not  given  by  Denton.  An  assumption  of  50  eV  (corre¬ 
sponding  to  mid  value  of  the  energy  range)  is  used  for  the  low  proton  temperature.  Once  the  data  is 
extracted  from  Denton  for  a  particular  turbulence  level  (decided  by  geomagnetic  Kp  index  value), 
an  interpolation  surface  is  fit  through  the  data  from  which  plasma  parameters  can  be  obtained  for 
any  orbit  local  time.  Low-charge  plasma  condition  is  simulated  using  Denton’s  data  for  Kp  —  1. 
For  simulating  high-charge  plasma  condition,  high-flux  ATS-6  data  [15],  given  in  Table  8,  have 
been  used. 


Table  8.  High- Charge  Plasma  Environment  for  GEO 


Particles 


Number  Density  (m  3)  Energy  (eV) 


electrons  2.36  x  105  16000 

protons  2.36  x  105  2.95  x  104 

Figures  14a  and  14b  show  voltage  and  equilibrium  currents,  respectively,  of  a  GEO  body  Under 
Low-charge  plasma  condition,  where  Plasma  Electron  Current  (PEC),  Plasma  Ion  Current  (PIC), 
Secondary  Electron  Current  (SEC),  Backscattered  Electron  Current  (BEC),  Photoelectric 
Current  (PC)  stand  for  plasma  electron  current,  respectively. 
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Unlike  LEO,  photoelectric  current  is  one  of  the  major  contributors  here,  and  its  absence  in 
shadow  region  results  in  large  dip  in  voltage.  The  body  potential  is  few  positive  volts  in 
sunlight.  On  the  other  hand,  for  high-charge  plasma  condition  described  by  Table  8,  body  voltage 
is  2.293  volts  under  sunlit  condition,  and  it  is  about  -50K  volts,  when  in  Earth’s  shadow. 


(a)  Body  Voltage  (b)  Logarithm  of  Absolute  Value  of 

Current 


Figure  14.  Charging  Characteristic  of  GEO  Group  1  (LAMR,  Low-Charge  Plasma 

Condition) 

Figure  15  shows  the  logarithm  of  magnitudes  of  perturbing  accelerations  for  GEO  orbits.  For  the 
LAMR  objects,  higher  harmonics  of  Earth’s  gravity  and  sun/moon  third  body  gravitational  pertur¬ 
bations  are  the  dominant  contributors  to  perturbing  acceleration;  for  the  HAMR  objects,  however, 
solar  radiation  pressure  also  becomes  one  of  the  major  contributors.  Lorentz  force  is  still  the  least 
dominating  factor  effecting  perturbing  acceleration.  The  sudden  jumps  in  Lorentz  perturbation  in 
Figure  15a  implies  that  the  object  is  in  Earth’s  shadow.  The  jumps  in  Lorentz  perturbation  in 
Figures  15b,  15c  are  because  of  magnetic  field  vector  nearly  aligning  up  with  the  relative  velocity 
vector. 
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(a)  GEO  Group  1 


(c)  GEO  Group  3 


(b)  GEO  Group  2 


. Sun  3rd  Body  Perturbation  LAMR 

■ . Moon  3rd  Body  Perturbation  LAMR 

. Solar  Radiation  Pressure  LAMR 

• . Lorentz  Force  LAMR 

Higher  Harmonics  of  Earth  Gravity  LAMR 
Sun  3rd  Body  Perturbation  HAMR 
Moon  3rd  Body  Perturbation  HAMR 

- Solar  Radiation  Pressure  HAMR 

- Lorentz  Force  HAMR 

Higher  Harmonics  of  Earth  Gravity  HAMR 


Figure  15.  Common  Logarithm  of  Perturbing  Acceleration  Magnitudes  for  GEO  Orbits, 

Low-Charge  Plasma  Condition 


Figures  16,  17,  18  show  in-track,  cross-track  and  radial  perturbations  of  Lorentz  force  included 
orbit  relative  to  Lorentz  force  excluded  orbit,  respectively  for  HAMR  objects.  Figures  16a,  16b 
show  the  in-track  perturbations  corresponding  to  low-charge  and  high-charge  plasma  conditions 
respectively.  Under  both  low-charge  and  high-charge  plasma  conditions,  lowly  inclined  group 
1  orbit  exhibits  secular  increasing  trend  in  in-track  direction,  whereas  higher  inclined  orbits  of 
group  2  and  group  3  exhibit  secular  decreasing  trend  in  in-track  plots.  In  addition  to  the  secular 
trend,  the  in-track  perturbations  also  exhibit  a  short  periodic  variation  with  a  period  of  one  day. 
The  increasing  secular  trend  in  group  1  in-track  perturbation  arises  because  of  decreasing  semi¬ 
major  axis  of  Lorentz  force  included  orbit  relative  to  Lorentz  force  excluded  orbit.  Similarly, 
the  decreasing  secular  trend  in  higher  inclined  in- track  plots  arise  out  of  increasing  semi-major 
axis  of  Lorentz  force  included  orbit  relative  to  Lorentz  force  excluded  orbit.  Since  the  in-track 
perturbations  for  group  1  orbit  are  positive,  object  in  Lorentz  force  included  orbit  lies  ahead  of 
object  in  Lorentz  force  excluded  orbit.  By  similar  logic,  object  in  higher  inclined  orbits  in  Lorentz 
force  included  orbits,  on  an  average,  lags  behind  object  in  Lorentz  force  excluded  orbits.  For 
group  1  orbit  Under  Low-charge  plasma  conditions,  the  order  of  in-track  perturbations  is  10  4  m 
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whereas  it  is  an  order  higher  for  group  2  and  group  3  orbits.  Under  High-charge  plasma  conditions, 
for  group  1  orbit,  the  order  of  in-track  perturbations  is  0.3  m,  which  is  three  orders  of  magnitude 
higher  than  group  2  and  group  3  orbits.  This  opposite  trend  for  group  1  orbit  under  High-charge 
plasma  condition  is  due  to  body  voltage  reaching  approximately  -50K  volts  in  shadow  Under  High- 
charge  plasma  condition;  as  the  simulations  are  run  near  equinox,  higher  inclined  orbits  do  not 
enter  Earth’s  shadow  region  and  hence  do  not  see  a  surge  in  Lorentz  force  values. 

Figures  17a,  17b  show  cross-track  perturbations  of  Lorentz  force  included  orbit  relative  to  Lorentz 
force  excluded  orbit  Under  Low-charge  and  high-charge  plasma  conditions  respectively.  The  cross¬ 
track  perturbation  plots  nearly  exhibit  a  zero-mean  periodic  variation  with  an  approximate  period 
of  one  day.  The  zero-mean  periodic  behavior  implies  that,  on  an  average,  there  is  no  change  in 
orbital  plane  due  to  Lorentz  force.  The  amplitudes  of  the  periodic  variations  in  the  cross-track 
perturbations  increase  with  time,  and  a  factor  responsible  for  this  increasing  behavior  is  the  de¬ 
creasing  RAAN  of  the  Lorentz  force  included  orbit  relative  to  Lorentz  force  excluded  orbit. 
Between  group  2  and  group  3  orbits,  group  3  orbits  exhibit  larger  amplitudes  of  variations  because 
of  higher  inclination.  For  low-charge  plasma  conditions,  cross-track  perturbations  are  of  the  order 
of  10  6  m  for  lowly  inclined  group  1  orbit,  whereas  it  is  two  orders  of  magnitude  higher  for  higher 
inclined  group  2  and  group  3  orbits.  For  group  1  orbit  under  high-charge  plasma  conditions,  cross¬ 
track  perturbations  are  of  the  order  of  1 0  3  m,  which  is  one  order  of  magnitude  higher  than  group  2 
and  group  3  orbits.  The  sharp  change  in  slope  of  group  1  orbit  under  high-charge  plasma  condition 
(near  .8  orbital  period  epoch)  corresponds  to  the  body  entering  Earth’s  shadow  region  for  the  first 
time. 

Figures  18a,  18b  show  radial  perturbations  of  Lorentz  force  included  orbit  relative  to  Lorentz 
force  excluded  orbit  Under  Low-charge  and  high-charge  plasma  conditions  respectively.  Group  1 
radial  perturbation  plots  clearly  exhibit  a  decreasing  secular  trend  superimposed  with  a  periodic 
variation.  The  secular  trend  arises  out  of  decreasing  semi-major  axis  of  Lorentz  force  included 
orbit  relative  to  Lorentz  force  excluded  orbit.  The  periodic  variation  in  radial  perturbations  of 
group  1  orbit  has  a  time-scale  of  approximately  one  day,  and  the  amplitude  of  variation  increases 
over  time.  Higher  inclined  group  2  and  group  3  orbits  have  a  small  secular  (increasing)  trend 
because  of  increasing  relative  semi-major  axis.  Group  2  and  group  3  orbits  also  have  a  periodic 
variation  of  one-day  period  with  increasing  amplitude.  Also,  group  3  orbits  have  larger  amplitudes 
of  variations  compared  to  group  2  orbits.  Under  Low-charge  plasma  conditions,  group  1  orbit  has 
radial  perturbations  of  the  order  of  10  6  m,  which  is  three  orders  of  magnitude  smaller  than  higher 
inclined  group  2  and  group  3  orbits.  On  the  contrary,  under  high-charge  plasma  conditions,  group 
1  orbit  has  perturbations  of  the  order  of  10  2  m,  which  is  an  order  of  magnitude  larger  than  higher 
inclined  groups. 
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(a)  In-Track  Perturbations,  Low-Charge 
Condition 


(b)  In-Track  Perturbations,  High-Charge 
Condition 


Figure  16.  In-Track  Perturbations  for  Lorentz  Force  Perturbed  GEO  Relative  to 
Lorentz  Force  Excluded  Orbit  (HAMR  Object) 


Number  of  Orbits 


(a)  Cross-track  Perturbations,  Low- 
charge  Condition 


(b)  Cross-track  Perturbations, 
High-charge  Condition 


Figure  17.  Cross-Track  Perturbations  for  Lorentz  Force  Perturbed  GEO  Relative  to 
Lorentz  Force  Excluded  Orbit  (HAMR  Object) 
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(a)  Radial  Perturbations,  Low-  (b)  Radial  Perturbations,  High- 

Charge  Condition  Charge  Condition 


Figure  18.  Radial  Perturbations  for  Lorentz  Force  Perturbed  GEO  Relative  to  Lorentz 

Force  Excluded  Orbit  (HAMR  Object) 


The  in-track,  cross-track  and  radial  perturbation  differences  appear  to  be  less  oscillatory  for  GEO 
when  compared  to  LEO.  The  primary  reason  for  this  difference  is  that  geosynchronous  body  com¬ 
pletes  only  four  orbits  whereas  low  Earth  body  completes  approximately  62  orbits  in  four  day 
simulation  period. 

It  is  interesting  to  investigate  the  orbital  element  differences  of  high-charge  group  1  orbit  be¬ 
cause  of  large  voltages  in  shadow.  Figure  19  shows  the  difference  in  orbital  elements  (between 
Lorentz  force  included  modeling  and  Lorentz  force  excluded  modeling)  for  group  1  orbits  un¬ 
der  high  charging  plasma  condition.  The  steep  changes  correspond  to  the  body  being  in  Earth’s 
shadow.  The  variations  in  difference  in  orbital  elements  (especially,  semi-major  axis  and  incli¬ 
nation)  in  between  the  jumps  is  small  relative  to  the  scale  of  the  jumps,  and  hence  the  apparent 
flattening. 


Number  of  Orbits 


Number  of  Orbits 


Number  of  Orbits 


Figure  19.  Difference  in  Semi-Major  Axis,  Eccentricity  and  Inclination  for  Cases  With 
and  Without  Lorentz  Force  Modeling.  GEO  Group  1  Under  High-Charge  Plasma 

Condition  (HAMR) 

Plots  corresponding  to  LAMR  objects  in  GEO  are  not  included  here  because  of  their  relative  small 
perturbation  values  when  compared  to  that  of  HAMR  objects.  Table  9,  however,  lists  the  orders  of 
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in-track,  cross-track  and  radial  perturbations  of  the  Lorentz  force  included  orbit  relative  to  Lorentz 
force  excluded  orbit  for  LAMR  objects. 


Table  9.  Orders  of  Magnitude  for  In-Track,  Cross-Track  and  Radial  Perturbations 
for  LAMR  Objects  in  GEO  Over  a  Period  of  Four  Days 


GEO  group 

In-track  perturbations  (m) 

Cross-track  perturbations 
(m) 

Radial  Perturbations  (m) 

Group  1 

1(T7 

ltr9 

1(T7 

Group  2 

1(T6 

1(T7 

10”7 

Group  3 

ltr6 

ltr7 

1(T6 

4.3  Comparison  of  Analytic  Expressions 

Table  10  compares  perturbation  values  obtained  using  analytic  expressions  developed  in  this  re¬ 
search  [Equations  (230)  -  (234)]  with  perturbation  values  obtained  using  analytic  expressions 
developed  by  Peng  (given  in  Appendix  C).  8  x  (x  =  a,  e,  i,  ft),  Q)  represents  difference  in  orbital 
elements  (after  one  orbital-period)  between  the  Lorentz-perturbed  and  Keplerian  orbits.  ‘A/  in  the 
tables  stands  for  analytic  expressions  developed  in  this  research  and  ‘Ap’  stands  for  analytic 
expressions  developed  by  Peng  et  al.  [14].  A  constant  voltage  of  -.4  volts  and  a  constant 
capacitance  of  7.72  x  10~9  F  are  assumed  for  the  LEO  object. 

Table  10.  Values  of  Orbital  Perturbations  for  LAMR  Objects  in  LEO  Under  Low- 

Charge  Plasma  Conditions 


Groups 

8a  (m) 

8e 

8i  (Deg.) 

8(0  (Deg.) 

8Q.  (Deg.) 

Group  1 

A,-:  -1.13E-9 
Ap.  -1.21E-9 

Ar:  1.02E-16 
Ap\  -1.26E-21 

Ar\  3.97E-11 
Ap\  4.27E-1 1 

A,-:  -3.95E-8 
Ap\  -3.98E-8 

A,-:  4.05E-8 

Ap\  4.08E-8 

Group  2 

A,-:  -6.58  E-7 
Ap\  -6.96E-7 

A,-:  5.72E-14 
Ap\  -1.37E-16 

Ar:  4.3  IE- 11 
Ap\  4.56E-11 

A,-:  -7.36E-10 
Ap\  1.75E-10 

Ar:  -4.82E-10 
Ap\  -4.83E-10 

Group  3 

A,-:  -6.55  E-7 
Ap.  -6.91E-7 

A,-:  5.67E-14 
Ap.  -1.38E-16 

Ar:  4.36E-11 
Ap\  4.60E-11 

Ar:  -9.30E-10 
Ap\  -3.14E-11 

Ar:  -4.97E-10 
Ap\  -4.98E-10 

As  can  be  seen  from  Table  10,  there  is  significant  difference  between  the  eccentricity  and  argument 
of  perigee  results.  Thorough  validation  is  needed  at  this  point. 

5.0  CONCLUSIONS 

A  detailed  modeling  of  plasma  electron  and  ion  currents,  backscattered  current,  secondary  current 
and  photoelectric  current  have  been  carried  out  for  computing  object  voltage  resulting  from  the 
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natural  space  environment.  An  approximate  space  plasma  environment  modeling  has  been 
achieved.  With  the  help  of  voltage  modeling  and  dipole  Earth  magnetic  field  modeling,  Lorentz 
forces  acting  as  orbital  perturbations  on  space  objects  has  been  computed.  The  Lorentz  forces 
are  induced  by  the  movement  of  the  objects  that  acquire  charge  in  the  natural  space  plasma 
environment.  The  derivations  have  been  made  under  the  assumption  of  a  spherical  conducting 
object.  Appropriate  models  for  Earth  gravitational  force,  Moon  and  Sun  third  body  gravitational 
forces,  solar  radiation  pressure  and  atmospheric  drag  have  also  been  included.  Sample 
simulations  have  been  made  to  illustrate  the  effect  of  the  Lorentz  force  perturbations 
representatively  for  high  area-to-mass  and  low  area-to-mass  ratio  objects,  that  are  spherical 
aluminum  conductors.  Simulations  in  geosyn-chronous  Earth  orbits  and  low  Earth  orbits  have 
been  made.  Both  low-charge  plasma  condition  and  high-charge  plasma  condition  have  been 
simulated. 

For  the  geosynchronous  region,  75°  longitude  position  has  been  selected  for  simulation.  Among 
low  Earth  orbits  (altitude  «  400  Km),  special  attention  has  been  paid  to  sun-synchronous  orbits. 
Effect  of  area-to-mass  ratio  and  inclination  have  been  studied  for  both  orbit  regimes.  Effect  of 
Earth’s  shadow  on  the  evolution  of  orbital  elements  has  been  studied.  Comparison  of  magnitudes 
of  different  perturbation  forces  has  been  done. 

For  low  Earth  low  area-to-mass  ratio  objects  under  low-charge  plasma  conditions,  in-track  per¬ 
turbations  exhibit  a  distinct  secular  increasing  trend,  cross-track  perturbations  exhibit  a  near  zero- 
mean  periodic  variation  with  increasing  amplitude,  and  the  radial  perturbations  of  highly  inclined 
orbits  (nearly  polar)  exhibit  a  secular  decreasing  trend.  For  low  Earth  low  area-to-mass  ratio  ob¬ 
jects  under  high-charge  plasma  conditions,  in-track  perturbations  exhibit  secular  increasing  trend 
for  group  2  (z,-  =  89°)  orbit,  and  secular  decreasing  trend  for  group  3  (sun-synchronous)  orbit;  the 
cross-track  and  radial  perturbations  exhibit  a  near- zero  mean  variation  with  increasing  amplitude 
values  .  For  geosynchronous  high  area-to-mass  ratio  objects,  under  both  low-charge  and  high- 
charge  plasma  conditions,  the  in-track  perturbations  exhibit  secular  decreasing  trend,  the  cross¬ 
track  perturbations  have  a  periodic  variation  with  near-zero  mean,  the  radial  perturbations  have 
a  distinct  secular  increasing  trend  for  lowly  inclined  (nearly  equatorial)  orbits  whereas  the  radial 
perturbations  have  a  small  secular  increasing  trend  for  higher  inclined  group  2  (z,  =  15°)  and  group 
3  (z'/  =  40°)  orbits. 

In  geosynchronous  region,  high-charge  plasma  condition  results  in  transition  of  voltage  from  few 
positive  volts  in  sunlight  to  approximately  -40  kY  in  shadow.  As  a  result  of  the  high  negative  volt¬ 
age,  in-track,  cross-track  and  radial  perturbations  for  lowly  inclined  high  area-to-mass  ratio  object 
increase  by  several  orders  of  magnitude  as  compared  to  low-charge  plasma  condition;  this  also 
brings  in  large  and  steep  changes  in  perturbations  of  orbital  elements.  In  low  Earth  region,  high- 
charge  plasma  condition  inside  auroral  oval  results  in  voltages  to  the  tune  of  -2.6  kV  in  sunlight 
and  to  the  tune  of  -60  kV  in  shadow.  This  results  in  in-track,  cross-track  and  radial  perturbations 
for  low  area-to-mass  ratio  object  increase  by  3-4  orders  of  magnitude  as  compared  to  low-charge 
plasma  condition. 

The  variational  equations  which  allow  for  the  derivation  of  the  secular  change  of  the  orbital  ele¬ 
ments  due  to  Lorentz  force  perturbations  have  been  significantly  extended.  A  thorough  validation 
and  in  depth  understanding  of  those  expressions  is  the  next  step  in  this  research.  The  analytical 
expressions  derived  from  the  variational  equations  allow  for  fast  long  term  propagation,  however 
further  research  is  needed. 
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APPENDIX  A  -  GENERAL  PERTURBATION  COEFFICIENTS 
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APPENDIX  B  -  LEO  PERTURBATION  COEFFICIENTS 
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APPENDIX  C  -  PERTURBATION  EXPRESSIONS  FROM  PENG  ET  AL.  see  [43] 
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APPENDIX  E  -  ABSTRACTS 


“Space  Debris  Charging  and  its  Effect  on  Orbit  Evolution”  published  in  AIAA/AAS 
Astrodynamics  Specialist  Conference:  With  the  increasing  number  of  debris  in  the  space 
environment  surrounding  Earth,  it  has  become  important  to  keep  track  of  the  orbits  of  these 
defunct  objects  so  as  to  avoid  collisions  with  active  satellites,  transiting  spacecrafts  or  other 
important  space  assets.  In  this  paper,  attention  has  been  paid  to  trajectory  evolution  of  debris 
in  low  Earth  orbit  and  geosynchronous  orbit  regions.  One  of  the  forces  effecting  the  trajec¬ 
tory  of  a  space  debris  is  the  Lorentz  force,  which  acts  when  a  charged  body  moves  through 
the  Earths  magnetosphere.  Because  of  continuous  bombardment  of  plasma  particles,  a  space 
debris  is  often  subject  to  charging.  Correct  modeling  of  Lorentz  force  requires  correct  mod¬ 
eling  of  magnetosphere  and  the  body  charge,  which  in  turn  depends  on  correct  modeling 
of  body  currents,  space-plasma  environment  and  body  capacitance.  This  research  involves 
modeling  of  in-space  charging  for  space  debris,  which  are  modeled  as  spherical  conductors. 
Simulations  incorporating  Lorentz  force  as  an  additional  perturbation  force  have  been  run  to 
evaluate  propagation  of  low  area-to-mass  ratio  and  high  area-to-mass  ratio  objects. 
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other  important  space  assets.  In  this  paper,  attention  has  been  paid  to  trajectory  evolution 
of  debris  in  low  Earth  orbit  and  geosynchronous  orbit  regions.  One  of  the  forces  effect¬ 
ing  the  trajectory  of  a  space  debris  is  the  Lorentz  force,  which  acts  when  a  charged  body 
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spherical  conductors.  Simulations  incorporating  Lorentz  force  as  an  additional  perturbation 
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LIST  OF  SYMBOLS,  ABBREVIATIONS,  AND  ACRONYMS 


Ae 

At 

B 

C 

Ibsc 

h 

h 

Iph 

Ise 

J 

L 

Te 

Ti 

V  sp 

z 

AMO 

AFRL 

BEC 

CGS 

CME 

DMSP 

DSCS 

GEO 

GOST 

HAMR 

HEOS 

IAGA 

IGRF 

IMP 

ISEE 

LAMR 

LEO 

LFM 

LVLH 

MHD 

MUSCAT 

NASA 


Electron  collection  area,  m2 
Ion  collection  area,  m 2 
Magnetic  Field,  T 
Capacitance,  F 
Backscattered  current,  A 
Plasma  electron  current,  A 
Plasma  ion  current,  A 
Photoelectric  current,  A 
Secondary  electron  current,  A 
Current  density,  A / m 2 
Debye  length,  m 
Electron  plasma  temperature,  K 
Ion  plasma  temperature,  K 
Object  speed,  m/s 
Atomic  number  of  object  material 
Air  Mass  Zero 

Air  Force  Research  Laboratory 
Backscatytered  Elecontron  Current 
Centimeter-  gram-  second 
Coronal  mass  ejection 
Defense  Meteorological  Satellite  Program 
Defense  Satellite  Communications  System 
Geosynchronous  Earth  orbit 
Gosudarstvennyy  Standart 
High  area-to-mass  ratio 
Highly  Eccentric  Orbit  Satellite 

International  Association  of  Geomagnetism  and  Aeronomy 

International  Geomagnetic  Reference  Field 

Interplanetary  Monitoring  Platform 

International  Sun-Earth  Explorer 

Low  area-to-mass  ratio 

Low  Earth  orbit 

Lyon-Fedder-Mobarry 

Local- vertical- local-horizontal 

Magnetohydrodynamic 

Multi-Utility  Spacecraft  Charging  Analysis  Tool 
National  Aeronautics  and  Space  Administration 


NASCAP-2K  NASA/ Air  Force  Spacecraft  Charging  Analysis  Program 
NASCAP-LEO  NASA  Charging  Analyzer  Program  for  Low-Earth  Orbit 
NASCAP-GEO  NASA  Charging  Analyzer  Program  for  Geosynchronous  Orbit 
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OGO 

Orbiting  Geophysical  Observatory 

PC 

Photoelectric  Current 

PEC 

Plasma  Electron  Current 

PIC 

Particle-in-cell/Plasma  Ion  Current 

POLAR 

Potentials  Of  Large  objects  in  the  Auroral  Region 

RAAN 

Right  ascension  of  ascending  node 

SCATHA 

Spacecraft  Charging  at  High  Altitudes 

SEC 

Secondary  Electron  Current 

SPARCS 

Spacecraft  Charging  Software 

SPIS 

Spacecraft  Plasma  Interaction  System 
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